Newer
Older
!-----------------------------------------------------------------
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source$ $Revision$
! MASDEV4_7 forcing 2006/10/16 14:23:23
!-----------------------------------------------------------------
! ###################
MODULE MODI_FORCING
! ###################
!
INTERFACE
!
SUBROUTINE FORCING ( PTSTEP, OUSERV, PRHODJ, PCORIOZ, &
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
PZHAT, PZZ, TPDTCUR, &
PUFRC_PAST, PVFRC_PAST, &
PUT, PVT, PWT, PTHT, PTKET, PRT, PSVT, &
PRUS, PRVS, PRWS, PRTHS, PRTKES, PRRS, PRSVS, &
KMI)
!
USE MODD_TIME, ONLY: DATE_TIME
!
REAL, INTENT(IN) :: PTSTEP ! time-step
LOGICAL , INTENT(IN) :: OUSERV ! Logical to use rv
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! ( rhod J ) = dry density
! for reference state * Jacobian of the GCS transformation.
REAL, DIMENSION(:,:), INTENT(IN) :: PCORIOZ ! f: Coriolis parameter
REAL, DIMENSION(:), INTENT(IN) :: PZHAT ! height level without orography
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height z
TYPE (DATE_TIME), INTENT(IN) :: TPDTCUR ! current date and time
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PUFRC_PAST, PVFRC_PAST
! ! forcing at previous time-step
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT,PVT,PWT,PTHT,PTKET
! wind, potential temperature and
! TKE at time t
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT ! moist variables at time t
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVT! scalar variables at time t
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS,PRVS,PRWS,PRTHS,PRTKES
! wind, potential temperature and
! TKE tendencies at time t
REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS ! moist variables at time t+1
REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRSVS! scalar variables at time t+1
!
INTEGER, INTENT(IN) :: KMI ! Model index
!
END SUBROUTINE FORCING
!
END INTERFACE
!
END MODULE MODI_FORCING
!
! ######################################################################
SUBROUTINE FORCING ( PTSTEP, OUSERV, PRHODJ, PCORIOZ, &
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
PZHAT, PZZ, TPDTCUR, &
PUFRC_PAST, PVFRC_PAST, &
PUT, PVT, PWT, PTHT, PTKET, PRT, PSVT, &
PRUS, PRVS, PRWS, PRTHS, PRTKES, PRRS, PRSVS, &
KMI)
! ######################################################################
!
!!*** *FORCING* - routine to compute the forced terms
!!
!! PURPOSE
!! -------
!! The routine prepares (linear interpolations) and integrates each
!! specified forcing terms. The forcing terms can be
!! - a geostrophic wind (u_frc, v_frc)
!! - a thermal wind (Dtheta/Dx, Dtheta/Dy)
!! - a tendency forcing (dth/dt, drv/dt )
!! - a vertical transport (w_frc)
!! - a newtonian relaxation (u_frc, v_frc, theta_frc, r_vfrc)
!!
!!** METHOD
!! ------
!! For its first call, the routine looks for a starting sounding with
!! a date_and_time immediately lower or close to that the current
!! date_and_time of the model. Then the temporal interpolation or extension
!! is performed according to the position of the current date_and_time
!! as compared to that of the forcing soundings. In case of non-flat
!! terrain, vertical interpolations are necessary because the forcing terms
!! are horizontally homogeneous.
!! All the necessary interpolations are linear. Integration of each forcing
!! term is enabled by a dedicated switch. The forced vertical motion is
!! computed for each prognostic variable with an upstream scheme applied to
!! the lagged fields. The forced advection of water vapor and the integration
!! of the thermal and geostrophic wind is time centered. If a forced
!! relaxation is enabled, a mask is defined to restrict the application of
!! the forcing.
!!
!! EXTERNAL
!! --------
!! Shuman functions (finite differences operators)
!! Upstream_z function (upstream selection of the vertical advective
!! tendency)
!! Temporal_lt function (compare 2 TYPEd date_and_time data)
!! Temporal_dist function (compute the number of seconds between
!! 2 TYPEd date_and_time data)
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!! Module MODD_CONF: declaration of configuration variables
!! LFLAT: tells if orography is present
!! L1D : tells if 1D model version is running
!! NVERB: level of printed information on the output listing
!! Module MODD_DYN: declaration of dynamic control variables
!! LCORIO: Coriolis parameter flag
!! Module MODD_FRC: declaration of the forcing variables
!! NFRC : number of forcing variables
!! TDTFRC: date of each forcing profile
!! XUFRC,XVFRC,XWFRC,XTHFRC,XRVFRC: large scale forcing variables
!! XGXTHFRC,XGYTHFRC: large scale gradient of Theta
!! XTENDTHFRC,XTENDRVFRC: large scale tendencies for Theta and Rv
!! Module MODD_LUNIT : contains logical unit names for all models
!! CLUOUT0 : name of output-listing
!! Module MODD_PARAMETERS: declaration of parameter variables
!! JPVEXT: define the number of marginal points out of the
!! physical domain along the vertical direction.
!! Module MODD_TIME: contains the structure of the TYPEd date_and_time
!!
!!
!! REFERENCE
!! ---------
!!
!! AUTHOR
!! ------
!! M. Georgelin * Laboratoire d'Aerologie*
!!
!! MODIFICATIONS
!! -------------
!! Original 13/12/95
!! 30/07/96 (Mari, Pinty, Suhre) restructuring
!! 18/11/96 J.-P. Pinty remove $n
!! 10/12/96 J.-P. Pinty reverse the loop order for vertical
!! interpolation in case of orography
!! 24/01/97 J.-P. Pinty add budget call
!! 24/01/98 P. Bechtold use tendency forcing instead of rv_advect
!! store large scale w
!! add SST and surface pressure forcing
!!
!! 01/2004 V. Masson surface externalization, removes SST
!! forcing
!! 06/2012 V. Masson Adds tendency of geostrophic wind itself to wind tendency
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODE_FM
USE MODE_IO_ll
!
USE MODD_CONF
USE MODD_DYN
USE MODD_FRC
USE MODD_LUNIT
USE MODD_PARAMETERS
USE MODD_TIME
USE MODD_BUDGET
!
USE MODI_SHUMAN
USE MODI_UPSTREAM_Z
USE MODI_TEMPORAL_LT
USE MODI_TEMPORAL_DIST
USE MODI_BUDGET
!
IMPLICIT NONE
!
!* 0.1 Declarations of dummy arguments :
!
REAL, INTENT(IN) :: PTSTEP ! time-step
LOGICAL , INTENT(IN) :: OUSERV ! Logical to use rv
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! ( rhod J ) = dry density
! for reference state * Jacobian of the GCS transformation.
REAL, DIMENSION(:,:), INTENT(IN) :: PCORIOZ ! f: Coriolis parameter
REAL, DIMENSION(:), INTENT(IN) :: PZHAT ! height level without orography
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height z
TYPE (DATE_TIME), INTENT(IN) :: TPDTCUR ! current date and time
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PUFRC_PAST, PVFRC_PAST
! ! forcing at previous time-step
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT,PVT,PWT,PTHT,PTKET
! wind, potential temperature and
! TKE at time t
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT ! moist variables at time t
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVT! scalar variables at time t
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS,PRVS,PRWS,PRTHS,PRTKES
! wind, potential temperature and
! TKE tendencies at time t
REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS ! moist variables at time t+1
REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRSVS! scalar variables at time t+1
!
INTEGER, INTENT(IN) :: KMI ! Model index
!
!* 0.2 Declarations of local variables
!
INTEGER :: IIU, IJU, IKU ! dimensions
INTEGER, SAVE :: JSX ! saved loop index
INTEGER :: JI, JJ, JK, JL, JXP! loop indexes
!
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZWF, ZUF, ZVF ! 3D forcing fields on
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTHF, ZRVF ! the model grid mesh
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZGXTHF, ZGYTHF ! at
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTENDTHF, ZTENDRVF ! time t
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDUF, ZDVF ! evolution of geostrophic wind
! ! during the time step
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZCOEF ! coefficient to take into
! ! account geostrophic wind evolution
! ! in wind tendencies
REAL,DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)) :: ZDUT, ZDVT ! variation of Wind components
! ! due to geostrophic wind evolution
!
!
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZZF, ZZA ! altitudes
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDZZ, ZRWCF ! dzz and ZWF contravar.
!
REAL, DIMENSION(SIZE(PUT,3)) :: ZXWFRC, ZXUFRC, ZXVFRC! 1D forcing fields
REAL, DIMENSION(SIZE(PUT,3)) :: ZXTHFRC, ZXRVFRC ! after
REAL, DIMENSION(SIZE(PUT,3)) :: ZXGXTHFRC, ZXGYTHFRC ! time
REAL, DIMENSION(SIZE(PUT,3)) :: ZXTENDTHFRC, ZXTENDRVFRC ! interpolation
REAL :: ZXPGROUNDFRC ! ground fields interpol.
!
LOGICAL, SAVE :: GSFIRSTCALL = .TRUE. ! control switch for the first call
!
REAL :: ZDZ, ZDT, ZALPHA ! height and time rate
REAL, SAVE :: ZSDTJX
REAL :: ZDZHAT_INV, ZDZHAT_INV_IKU ! Inverse of the vertical mesh size
!
INTEGER :: ILUOUT0 ! Logical unit number for output-listing
INTEGER :: IRESP ! Return code of FM-routines
!
LOGICAL,DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)) :: GRELAX_MASK_FRC
!
!----------------------------------------------------------------------------
!
IIU=SIZE(PUT,1)
IJU=SIZE(PUT,2)
IKU=SIZE(PUT,3)
!
!* 1. PREPARATION OF FORCING
! ----------------------
!
IF (GSFIRSTCALL) THEN
!
GSFIRSTCALL = .FALSE.
!
CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT0,IRESP)
!
!* 1.1 printout number of forcing profiles
!
WRITE(UNIT=ILUOUT0,FMT='(" THERE ARE ",I2," FORCING SOUNDINGS AT:")') NFRC
DO JSX = 1 , NFRC
WRITE(UNIT=ILUOUT0,FMT='(F9.0, "s, date:", I3, "/", I3, "/", I5)') &
TDTFRC(JSX)%TIME, &
TDTFRC(JSX)%TDATE%DAY, &
TDTFRC(JSX)%TDATE%MONTH, &
TDTFRC(JSX)%TDATE%YEAR
END DO
!
!* 1.2 find first sounding to be used
!
JSX = 0
IF( TEMPORAL_LT ( TPDTCUR, TDTFRC(1) ) ) THEN
WRITE(UNIT=ILUOUT0,FMT='(" THE INITIAL FORCING FIELDS ARE NULL ")')
ELSE IF( .NOT. TEMPORAL_LT ( TPDTCUR, TDTFRC(NFRC) ) ) THEN
WRITE(UNIT=ILUOUT0,FMT='(" THE FORCING FIELDS WILL REMAIN STATIONARY ")')
ELSE
!
TIM_FOR: DO JI = NFRC-1, 1, -1
JSX = JI
IF( .NOT. TEMPORAL_LT ( TPDTCUR, TDTFRC(JSX) ) ) EXIT TIM_FOR
END DO TIM_FOR
!
WRITE(UNIT=ILUOUT0,FMT='(" THE INITIAL FORCING FIELDS ARE INTERPOLATED" , &
& " IN TIME STARTING FROM THE SOUNDING NUMBER ",I2)') JSX
JSX = JSX - 1
END IF
!
!* 1.3 printout all forcing field
!
IF (NVERB >= 5) THEN
!
WRITE(UNIT=ILUOUT0,FMT='(A)') &
"FORCING: the following forcing fields will be used:"
!
WRITE(UNIT=ILUOUT0,FMT='(A)') &
"XUFRC: geostrophic wind in X"
DO JK = 1, IKU
WRITE(UNIT=ILUOUT0,FMT='(I10,99(/8F10.2))') &
JK, (XUFRC(JK,JL), JL=1, NFRC)
END DO
!
WRITE(UNIT=ILUOUT0,FMT='(A)') &
"XVFRC: geostrophic wind in Y"
DO JK = 1, IKU
WRITE(UNIT=ILUOUT0,FMT='(I10,99(/8F10.2))') &
JK, (XVFRC(JK,JL), JL=1, NFRC)
END DO
!
WRITE(UNIT=ILUOUT0,FMT='(A)') &
"XTHFRC: THETA for relaxation"
DO JK = 1, IKU
WRITE(UNIT=ILUOUT0,FMT='(I10,99(/8F10.2))') &
JK, (XTHFRC(JK,JL), JL=1, NFRC)
END DO
!
WRITE(UNIT=ILUOUT0,FMT='(A)') &
"XRVFRC: RV for relaxation"
DO JK = 1, IKU
WRITE(UNIT=ILUOUT0,FMT='(I10,99(/8F10.6))') &
JK, (XRVFRC(JK,JL), JL=1, NFRC)
END DO
!
WRITE(UNIT=ILUOUT0,FMT='(A)') &
"XWFRC: vertical transport velocity"
DO JK = 1, IKU
WRITE(UNIT=ILUOUT0,FMT='(I10,99(/8F10.7))') &
JK, (XWFRC(JK,JL), JL=1, NFRC)
END DO
!
WRITE(UNIT=ILUOUT0,FMT='(A)') &
"XGXTHFRC: thermal wind advection in X"
DO JK = 1, IKU
WRITE(UNIT=ILUOUT0,FMT='(I10,99(/8E10.3))') &
JK, (XGXTHFRC(JK,JL), JL=1, NFRC)
END DO
!
WRITE(UNIT=ILUOUT0,FMT='(A)') &
"XGYTHFRC: thermal wind advection in Y"
DO JK = 1, IKU
WRITE(UNIT=ILUOUT0,FMT='(I10,99(/8E10.3))') &
JK, (XGYTHFRC(JK,JL), JL=1, NFRC)
END DO
!
WRITE(UNIT=ILUOUT0,FMT='(A)') &
"XTENDTHFRC: Theta tendency"
DO JK = 1, IKU
WRITE(UNIT=ILUOUT0,FMT='(I10,99(/8E10.3))') &
JK, (XTENDTHFRC(JK,JL), JL=1, NFRC)
END DO
!
WRITE(UNIT=ILUOUT0,FMT='(A)') &
"XTENDRVFRC: humidity tendency"
DO JK = 1, IKU
WRITE(UNIT=ILUOUT0,FMT='(I10,99(/8E10.3))') &
JK, (XTENDRVFRC(JK,JL), JL=1, NFRC)
END DO
!
WRITE(UNIT=ILUOUT0,FMT='(A)') &
"XPGROUNDFRC: SURF PRESSURE FORCING"
WRITE(UNIT=ILUOUT0,FMT='(10X,99(/8E10.3))') &
(XPGROUNDFRC(JL), JL=1, NFRC)
!
END IF
!
END IF
!
!* 2. INTERPOLATE IN TIME
! -------------------
!
IF( TEMPORAL_LT ( TPDTCUR, TDTFRC(1) ) ) THEN
ZXUFRC(:) = XUFRC(:,1)
ZXVFRC(:) = XVFRC(:,1)
ZXWFRC(:) = XWFRC(:,1)
ZXTHFRC(:) = XTHFRC(:,1)
ZXRVFRC(:) = XRVFRC(:,1)
ZXTENDTHFRC(:) = XTENDTHFRC(:,1)
ZXTENDRVFRC(:) = XTENDRVFRC(:,1)
ZXGXTHFRC(:) = XGXTHFRC(:,1)
ZXGYTHFRC(:) = XGYTHFRC(:,1)
ZXPGROUNDFRC = XPGROUNDFRC(1)
ELSE IF ( .NOT. TEMPORAL_LT ( TPDTCUR, TDTFRC(NFRC) ) ) THEN
ZXUFRC(:) = XUFRC(:,NFRC)
ZXVFRC(:) = XVFRC(:,NFRC)
ZXWFRC(:) = XWFRC(:,NFRC)
ZXTHFRC(:) = XTHFRC(:,NFRC)
ZXRVFRC(:) = XRVFRC(:,NFRC)
ZXTENDTHFRC(:) = XTENDTHFRC(:,NFRC)
ZXTENDRVFRC(:) = XTENDRVFRC(:,NFRC)
ZXGXTHFRC(:) = XGXTHFRC(:,NFRC)
ZXGYTHFRC(:) = XGYTHFRC(:,NFRC)
ZXPGROUNDFRC = XPGROUNDFRC(NFRC)
ELSE
JXP = JSX + 1
IF( .NOT. TEMPORAL_LT ( TPDTCUR, TDTFRC(JXP) ) ) THEN
JSX = JSX +1
JXP= JSX +1
CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT0,IRESP)
WRITE(UNIT=ILUOUT0,FMT='(" THE FORCING FIELDS ARE INTERPOLATED NOW" ,&
& " BETWEEN SOUNDING NUMBER ",I2," AND SOUNDING NUMBER ",I2)') JSX,JXP
CALL TEMPORAL_DIST ( TDTFRC(JXP)%TDATE%YEAR,TDTFRC(JXP)%TDATE%MONTH, &
TDTFRC(JXP)%TDATE%DAY ,TDTFRC(JXP)%TIME, &
TDTFRC(JSX)%TDATE%YEAR ,TDTFRC(JSX)%TDATE%MONTH, &
TDTFRC(JSX)%TDATE%DAY ,TDTFRC(JSX)%TIME, &
ZSDTJX )
END IF
!
CALL TEMPORAL_DIST ( TPDTCUR%TDATE%YEAR ,TPDTCUR%TDATE%MONTH, &
TPDTCUR%TDATE%DAY ,TPDTCUR%TIME, &
TDTFRC(JSX)%TDATE%YEAR,TDTFRC(JSX)%TDATE%MONTH, &
TDTFRC(JSX)%TDATE%DAY ,TDTFRC(JSX)%TIME, &
ZDT )
!
ZALPHA = ZDT / ZSDTJX
!
ZXUFRC(:) = XUFRC(:,JSX) +(XUFRC(:,JXP)-XUFRC(:,JSX))*ZALPHA
ZXVFRC(:) = XVFRC(:,JSX) +(XVFRC(:,JXP)-XVFRC(:,JSX))*ZALPHA
ZXWFRC(:) = XWFRC(:,JSX) +(XWFRC(:,JXP)-XWFRC(:,JSX))*ZALPHA
ZXTHFRC(:) = XTHFRC(:,JSX) +(XTHFRC(:,JXP)-XTHFRC(:,JSX))*ZALPHA
ZXRVFRC(:) = XRVFRC(:,JSX) +(XRVFRC(:,JXP)-XRVFRC(:,JSX))*ZALPHA
ZXTENDTHFRC(:) = XTENDTHFRC(:,JSX)+(XTENDTHFRC(:,JXP)-XTENDTHFRC(:,JSX))*ZALPHA
ZXTENDRVFRC(:) = XTENDRVFRC(:,JSX)+(XTENDRVFRC(:,JXP)-XTENDRVFRC(:,JSX))*ZALPHA
ZXGXTHFRC(:) = XGXTHFRC(:,JSX)+(XGXTHFRC(:,JXP)-XGXTHFRC(:,JSX))*ZALPHA
ZXGYTHFRC(:) = XGYTHFRC(:,JSX)+(XGYTHFRC(:,JXP)-XGYTHFRC(:,JSX))*ZALPHA
ZXPGROUNDFRC = XPGROUNDFRC(JSX) +(XPGROUNDFRC(JXP)-XPGROUNDFRC(JSX))*ZALPHA
!
END IF
!
!
!* 3. INTERPOLATE IN SPACE
! --------------------
!
ALLOCATE(ZUF(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)))
ALLOCATE(ZVF(SIZE(PVT,1),SIZE(PVT,2),SIZE(PVT,3)))
ALLOCATE(ZWF(SIZE(PWT,1),SIZE(PWT,2),SIZE(PWT,3)))
ALLOCATE(ZTHF(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)))
ALLOCATE(ZRVF(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)))
ALLOCATE(ZGXTHF(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)))
ALLOCATE(ZGYTHF(SIZE(PVT,1),SIZE(PVT,2),SIZE(PVT,3)))
ALLOCATE(ZTENDTHF(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)))
ALLOCATE(ZTENDRVF(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)))
ALLOCATE(ZDUF(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)))
ALLOCATE(ZDVF(SIZE(PVT,1),SIZE(PVT,2),SIZE(PVT,3)))
!
IF (LFLAT) THEN
!
!* 3.1 flat terrain case
!
ZUF(:,:,:) = SPREAD( SPREAD( ZXUFRC(:),1,IIU ) ,2,IJU )
ZVF(:,:,:) = SPREAD( SPREAD( ZXVFRC(:),1,IIU ) ,2,IJU )
ZWF(:,:,:) = SPREAD( SPREAD( ZXWFRC(:),1,IIU ) ,2,IJU )
ZTHF(:,:,:) = SPREAD( SPREAD( ZXTHFRC(:),1,IIU ) ,2,IJU )
ZRVF(:,:,:) = SPREAD( SPREAD( ZXRVFRC(:),1,IIU ) ,2,IJU )
ZTENDTHF(:,:,:) = SPREAD( SPREAD( ZXTENDTHFRC(:),1,IIU ),2,IJU )
ZTENDRVF(:,:,:) = SPREAD( SPREAD( ZXTENDRVFRC(:),1,IIU ),2,IJU )
ZGXTHF(:,:,:) = SPREAD( SPREAD( ZXGXTHFRC(:),1,IIU ),2,IJU )
ZGYTHF(:,:,:) = SPREAD( SPREAD( ZXGYTHFRC(:),1,IIU ),2,IJU )
ELSE
!
!* 3.2 case with orography
!
ALLOCATE(ZZA(SIZE(PWT,1),SIZE(PWT,2),SIZE(PWT,3)))
ALLOCATE(ZZF(SIZE(PWT,1),SIZE(PWT,2),SIZE(PWT,3)))
ZZA(:,:,:) = MZF(1,IKU,1, PZZ(:,:,:) )
ZZA(:,:,IKU) = 2.0*PZZ(:,:,IKU) - ZZA(:,:,IKU-1)
ZDZHAT_INV_IKU = 1.0 / ( PZHAT(IKU)-PZHAT(IKU-1) )
!
ZZF(:,:,:) = MXM( ZZA(:,:,:) )
ZZF(1,:,:) = 2.0*ZZA(1,:,:) - ZZF(2,:,:)
!
DO JL=1,IKU-1
ZDZHAT_INV = 1.0 / ( PZHAT(JL+1)-PZHAT(JL) )
DO JK=1,IKU
DO JJ=1,IJU
DO JI=1,IIU
IF( ZZF(JI,JJ,JK) >= PZHAT(JL).AND.ZZF(JI,JJ,JK) <= PZHAT(JL+1) ) THEN
ZDZ = (ZZF(JI,JJ,JK)-PZHAT(JL)) * ZDZHAT_INV
ZUF(JI,JJ,JK) = ZXUFRC(JL+1)*ZDZ + ZXUFRC(JL)*(1-ZDZ)
ELSE IF( ZZF(JI,JJ,JK) > PZHAT(IKU) ) THEN
ZDZ = (ZZF(JI,JJ,JK)-PZHAT(IKU)) * ZDZHAT_INV_IKU
ZUF(JI,JJ,JK) = ZXUFRC(IKU)*ZDZ + ZXUFRC(IKU-1)*(1-ZDZ)
END IF
END DO
END DO
END DO
END DO
!
ZZF(:,:,:) = MYM( ZZA(:,:,:) )
ZZF(:,1,:) = 2.0*ZZA(:,1,:)-ZZF(:,2,:)
!
DO JL=1,IKU-1
ZDZHAT_INV = 1.0 / ( PZHAT(JL+1)-PZHAT(JL) )
DO JK=1,IKU
DO JJ=1,IJU
DO JI=1,IIU
IF( ZZF(JI,JJ,JK) >= PZHAT(JL).AND.ZZF(JI,JJ,JK) <= PZHAT(JL+1) ) THEN
ZDZ = (ZZF(JI,JJ,JK)-PZHAT(JL)) * ZDZHAT_INV
ZVF(JI,JJ,JK) = ZXVFRC(JL+1)*ZDZ + ZXVFRC(JL)*(1-ZDZ)
ELSE IF( ZZF(JI,JJ,JK) > PZHAT(IKU) ) THEN
ZDZ = (ZZF(JI,JJ,JK)-PZHAT(IKU)) * ZDZHAT_INV_IKU
ZVF(JI,JJ,JK) = ZXVFRC(IKU)*ZDZ + ZXVFRC(IKU-1)*(1-ZDZ)
END IF
END DO
END DO
END DO
END DO
DO JL=1,IKU-1
ZDZHAT_INV = 1.0 / ( PZHAT(JL+1)-PZHAT(JL) )
DO JK=1,IKU
DO JJ=1,IJU
DO JI=1,IIU
IF( ZZF(JI,JJ,JK) >= PZHAT(JL).AND.ZZF(JI,JJ,JK) <= PZHAT(JL+1) ) THEN
ZDZ = (ZZF(JI,JJ,JK)-PZHAT(JL)) * ZDZHAT_INV
ZWF(JI,JJ,JK) = ZXWFRC(JL+1)*ZDZ + ZXWFRC(JL)*(1-ZDZ)
ELSE IF( ZZF(JI,JJ,JK) > PZHAT(IKU) ) THEN
ZDZ = (ZZF(JI,JJ,JK)-PZHAT(IKU)) * ZDZHAT_INV_IKU
ZWF(JI,JJ,JK) = ZXWFRC(IKU)*ZDZ + ZXWFRC(IKU-1)*(1-ZDZ)
END IF
END DO
END DO
END DO
END DO
!
ZZF(:,:,:) = MZF(1,IKU,1, PZZ(:,:,:) )
ZZF(:,:,IKU) = 2.0*PZZ(:,:,IKU)-ZZF(:,:,IKU-1)
!
DO JL=1,IKU-1
ZDZHAT_INV = 1.0 / ( PZHAT(JL+1)-PZHAT(JL) )
DO JK=1,IKU
DO JJ=1,IJU
DO JI=1,IIU
IF( ZZF(JI,JJ,JK) >= PZHAT(JL).AND.ZZF(JI,JJ,JK) <= PZHAT(JL+1) ) THEN
ZDZ = (ZZF(JI,JJ,JK)-PZHAT(JL)) * ZDZHAT_INV
ZTHF(JI,JJ,JK) = ZXTHFRC(JL+1)*ZDZ + ZXTHFRC(JL)*(1-ZDZ)
ZRVF(JI,JJ,JK) = ZXRVFRC(JL+1)*ZDZ + ZXRVFRC(JL)*(1-ZDZ)
ZGXTHF(JI,JJ,JK) = ZXGXTHFRC(JL+1)*ZDZ + ZXGXTHFRC(JL)*(1-ZDZ)
ZGYTHF(JI,JJ,JK) = ZXGYTHFRC(JL+1)*ZDZ + ZXGYTHFRC(JL)*(1-ZDZ)
ZTENDTHF(JI,JJ,JK) = ZXTENDTHFRC(JL+1)*ZDZ + ZXTENDTHFRC(JL)*(1-ZDZ)
ZTENDRVF(JI,JJ,JK) = ZXTENDRVFRC(JL+1)*ZDZ + ZXTENDRVFRC(JL)*(1-ZDZ)
ELSE IF( ZZF(JI,JJ,JK) > PZHAT(IKU) ) THEN
ZDZ = (ZZF(JI,JJ,JK)-PZHAT(IKU)) * ZDZHAT_INV_IKU
ZTHF(JI,JJ,JK) = ZXTHFRC(IKU)*ZDZ + ZXTHFRC(IKU-1)*(1-ZDZ)
ZRVF(JI,JJ,JK) = ZXRVFRC(IKU)*ZDZ + ZXRVFRC(IKU-1)*(1-ZDZ)
ZGXTHF(JI,JJ,JK) = ZXGXTHFRC(IKU)*ZDZ + ZXGXTHFRC(IKU-1)*(1-ZDZ)
ZGYTHF(JI,JJ,JK) = ZXGYTHFRC(IKU)*ZDZ + ZXGYTHFRC(IKU-1)*(1-ZDZ)
ZTENDTHF(JI,JJ,JK) = ZXTENDTHFRC(IKU)*ZDZ + ZXTENDTHFRC(IKU-1)*(1-ZDZ)
ZTENDRVF(JI,JJ,JK) = ZXTENDRVFRC(IKU)*ZDZ + ZXTENDRVFRC(IKU-1)*(1-ZDZ)
END IF
END DO
END DO
END DO
END DO
END IF
!
!
! under the ground, forcings do not exist.
!
DO JK=1,JPVEXT
ZUF(:,:,JK) = 0.
ZVF(:,:,JK) = 0.
ZWF(:,:,JK) = 0.
ZTHF(:,:,JK) = 0.
ZRVF(:,:,JK) = 0.
ZGXTHF(:,:,JK) = 0.
ZGYTHF(:,:,JK) = 0.
ZTENDTHF(:,:,JK) = 0.
ZTENDRVF(:,:,JK) = 0.
END DO
!
! store large scale w in module to be used later
! in convection scheme
XWTFRC(:,:,:) = ZWF(:,:,:)
!
!* computes evolution of forcing wind
WHERE(PUFRC_PAST==XUNDEF) PUFRC_PAST = ZUF(:,:,:)
WHERE(PVFRC_PAST==XUNDEF) PVFRC_PAST = ZVF(:,:,:)
!
ZDUF(:,:,:) = ZUF(:,:,:) - PUFRC_PAST(:,:,:)
ZDVF(:,:,:) = ZVF(:,:,:) - PVFRC_PAST(:,:,:)
!
!
IF (.NOT.LFLAT) THEN
!
DEALLOCATE(ZZA)
DEALLOCATE(ZZF)
!
END IF
!
!
!* 4. INTEGRATION OF THE FORCINGS IN THE SOURCES
! ------------------------------------------
!
ALLOCATE(ZDZZ(SIZE(PWT,1),SIZE(PWT,2),SIZE(PWT,3)))
ALLOCATE(ZRWCF(SIZE(PWT,1),SIZE(PWT,2),SIZE(PWT,3)))
!
!* 4.1 integration of vertical motion (upstream scheme)
!
IF (LVERT_MOTION_FRC) THEN
ZDZZ(:,:,:) = DZM(1,IKU,1,MZF(1,IKU,1,PZZ(:,:,:)))
ZDZZ(:,:,IKU) = PZZ(:,:,IKU) - PZZ(:,:,IKU-1) ! same delta z in IKU and IKU -1
!
ZRWCF(:,:,:) = ZWF(:,:,:) * MZM(1,IKU,1,PRHODJ(:,:,:)) / ZDZZ(:,:,:)
ZRWCF(:,:,1) = - ZRWCF(:,:,3) ! Mirror hypothesis
!
! forced vertical transport of U and V
!
ZDZZ(:,:,:) = MXF(ZRWCF(:,:,:)) *DZM(1,IKU,1,PUT(:,:,:))
PRUS(:,:,:) = PRUS(:,:,:) - UPSTREAM_Z(ZDZZ(:,:,:),ZRWCF(:,:,:))
ZDZZ(:,:,:) = MYF(ZRWCF(:,:,:)) *DZM(1,IKU,1,PVT(:,:,:))
PRVS(:,:,:) = PRVS(:,:,:) - UPSTREAM_Z(ZDZZ(:,:,:),ZRWCF(:,:,:))
!
! forced vertical transport of W
!
IF( .NOT.L1D ) THEN
ZDZZ(:,:,:) = MZF(1,IKU,1,ZRWCF(:,:,:)) *DZF(1,IKU,1,PWT(:,:,:))
PRWS(:,:,:) = PRWS(:,:,:) - UPSTREAM_Z(ZDZZ(:,:,:),ZRWCF(:,:,:))
END IF
!
! forced vertical transport of THETA
!
ZDZZ(:,:,:) = ZRWCF(:,:,:) *DZM(1,IKU,1,PTHT(:,:,:))
PRTHS(:,:,:) = PRTHS(:,:,:) - UPSTREAM_Z(ZDZZ(:,:,:),ZRWCF(:,:,:))
!
! forced vertical transport of TKE (if allocated)
!
IF( SIZE(PTKET) == SIZE(ZDZZ) ) THEN
ZDZZ(:,:,:) = ZRWCF(:,:,:) *DZM(1,IKU,1,PTKET(:,:,:))
PRTKES(:,:,:) = PRTKES(:,:,:) - UPSTREAM_Z(ZDZZ(:,:,:),ZRWCF(:,:,:))
END IF
!
! forced vertical transport of water variables
!
DO JL = 1 , SIZE(PRRS,4)
ZDZZ(:,:,:) = ZRWCF(:,:,:) *DZM(1,IKU,1,PRT(:,:,:,JL))
PRRS(:,:,:,JL) = PRRS(:,:,:,JL) - UPSTREAM_Z(ZDZZ(:,:,:),ZRWCF(:,:,:))
END DO
!
! forced vertical transport of scalar variables
!
DO JL = 1 , SIZE(PRSVS,4)
ZDZZ(:,:,:) = ZRWCF(:,:,:) *DZM(1,IKU,1,PSVT(:,:,:,JL))
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
PRSVS(:,:,:,JL) = PRSVS(:,:,:,JL) - UPSTREAM_Z(ZDZZ(:,:,:),ZRWCF(:,:,:))
END DO
!
END IF
!
!* 4.2 integration of the tendency forcing for th and rv
!
IF ( LTEND_THRV_FRC ) THEN
PRTHS(:,:,:) = PRTHS(:,:,:) + PRHODJ(:,:,:) * ZTENDTHF(:,:,:)
IF ( OUSERV ) THEN
PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRHODJ(:,:,:) * ZTENDRVF(:,:,:)
END IF
END IF
!
!* 4.3 integration of the thermal and geostrophic wind
!
IF( LCORIO ) THEN
!
! thermal wind advection
!
IF (LGEOST_TH_FRC) THEN
PRTHS(:,:,:) = PRTHS(:,:,:) - PRHODJ(:,:,:)*(MXF(PUT(:,:,:))*ZGXTHF(:,:,:) &
+ MYF(PVT(:,:,:))*ZGYTHF(:,:,:) )
END IF
!
! geostrophic wind for U and V due to large scale pressure gradients
!
IF (LGEOST_UV_FRC) THEN
! Adds pressure force (in the form of geostrophic wind) to U component
PRUS(:,:,:) = PRUS(:,:,:) &
- MXM( MYF(ZVF(:,:,:))*PRHODJ(:,:,:)*SPREAD(PCORIOZ(:,:),3,IKU))
! Adds pressure force (in the form of geostrophic wind) to V component
PRVS(:,:,:) = PRVS(:,:,:) &
+ MYM( MXF(ZUF(:,:,:))*PRHODJ(:,:,:)*SPREAD(PCORIOZ(:,:),3,IKU))
! adds tendency of geostrophic wind to force wind in the free troposphere to
! follow the geostrophic wind when the latter changes.
! When winds differs from the geotrophic wind, the impact of this tendency is reduced.
ALLOCATE(ZCOEF(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)))
ZCOEF(:,:,:) = (MXF(PUT **2)+MYF(PVT **2)) &
/MAX(MXF(PUFRC_PAST**2)+MYF(PVFRC_PAST**2), 1.E-3)
!
ZCOEF(:,:,:) = MIN(1.,SQRT(ZCOEF))
!
ZDUT(:,:,:) = ZDUF(:,:,:) * MXM(ZCOEF)
ZDVT(:,:,:) = ZDVF(:,:,:) * MYM(ZCOEF)
PRUS(:,:,:) = PRUS(:,:,:) + ZDUT(:,:,:) * MXM(PRHODJ) / PTSTEP
PRVS(:,:,:) = PRVS(:,:,:) + ZDVT(:,:,:) * MYM(PRHODJ) / PTSTEP
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
!
!
! Takes into acount the Coriolis force due to this evolution
PRUS(:,:,:) = PRUS(:,:,:) &
+ MXM( MYF(ZDVT(:,:,:))*PRHODJ(:,:,:)*SPREAD(PCORIOZ(:,:),3,IKU))
PRVS(:,:,:) = PRVS(:,:,:) &
- MYM( MXF(ZDUT(:,:,:))*PRHODJ(:,:,:)*SPREAD(PCORIOZ(:,:),3,IKU))
!
DEALLOCATE(ZCOEF)
END IF
!
END IF
!
! stores new forcing wind
PUFRC_PAST(:,:,:) = ZUF(:,:,:)
PVFRC_PAST(:,:,:) = ZVF(:,:,:)
!
!
!* 4.4 integration of the thermal, moisture and wind relaxation
!
IF( LRELAX_THRV_FRC .OR. LRELAX_UV_FRC ) THEN
!
ZDZZ(:,:,:) = DZM(1,IKU,1,MZF(1,IKU,1,PZZ(:,:,:)))
ZDZZ(:,:,IKU) = PZZ(:,:,IKU) - PZZ(:,:,IKU-1)
!
! define the mask where the relaxation is to be applied
!
SELECT CASE (CRELAX_HEIGHT_TYPE)
CASE ("FIXE")
GRELAX_MASK_FRC(:,:,:) = .TRUE.
CASE ("THGR")
CALL DEFINE_RELAX_FORCING(GRELAX_MASK_FRC,PTHT,ZDZZ,JPHEXT,JPVEXT)
CASE DEFAULT
! the following error should not occur, since tests are made earlier
!callabortstop
CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
CALL ABORT
STOP "Error in FORCING: wrong CRELAX_HEIGHT_TYPE option."
END SELECT
WHERE ( MZF(1,IKU,1,PZZ(:,:,:)) .LE. XRELAX_HEIGHT_FRC )
GRELAX_MASK_FRC = .FALSE.
END WHERE
!
IF( LRELAX_THRV_FRC ) THEN
!
! apply THETA relaxation
!
WHERE( GRELAX_MASK_FRC )
PRTHS(:,:,:) = PRTHS(:,:,:) - PRHODJ(:,:,:)*(PTHT(:,:,:)-ZTHF(:,:,:)) &
/ XRELAX_TIME_FRC
END WHERE
!
! apply humidity relaxation
!
IF( OUSERV ) THEN
WHERE( GRELAX_MASK_FRC )
PRRS(:,:,:,1) = PRRS(:,:,:,1) &
/ XRELAX_TIME_FRC
END WHERE
!
END IF
!
END IF
!
IF( LRELAX_UV_FRC ) THEN
!
! apply UV relaxation
!
WHERE( GRELAX_MASK_FRC )
PRUS(:,:,:) = PRUS(:,:,:) - MXM(PRHODJ(:,:,:))*(PUT(:,:,:)-ZUF(:,:,:)) &
PRVS(:,:,:) = PRVS(:,:,:) - MYM(PRHODJ(:,:,:))*(PVT(:,:,:)-ZVF(:,:,:)) &
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
/ XRELAX_TIME_FRC
END WHERE
!
END IF
!
END IF
!
!
!* 4.6 ground pressure forcing
!
IF( LPGROUND_FRC ) THEN
!
! to be implemented as a function of ZXPGROUNFRC
!
END IF
!
!
!* 5. BUDGET CALLS
! ------------
!
!
IF (LBUDGET_U) CALL BUDGET (PRUS,1,'FRC_BU_RU')
IF (LBUDGET_V) CALL BUDGET (PRVS,2,'FRC_BU_RV')
IF (LBUDGET_W) CALL BUDGET (PRWS,3,'FRC_BU_RW')
IF (LBUDGET_TH) CALL BUDGET (PRTHS,4,'FRC_BU_RTH')
IF (LBUDGET_TKE) CALL BUDGET (PRTKES,5,'FRC_BU_RTKE')
IF (LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1),6,'FRC_BU_RRV')
IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2),7,'FRC_BU_RRC')
IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:,3),8,'FRC_BU_RRR')
IF (LBUDGET_RI) CALL BUDGET (PRRS(:,:,:,4),9,'FRC_BU_RRI')
IF (LBUDGET_RS) CALL BUDGET (PRRS(:,:,:,5),10,'FRC_BU_RRS')
IF (LBUDGET_RG) CALL BUDGET (PRRS(:,:,:,6),11,'FRC_BU_RRG')
IF (LBUDGET_RH) CALL BUDGET (PRRS(:,:,:,7),12,'FRC_BU_RRH')
IF (LBUDGET_SV) THEN
DO JL = 1 , SIZE(PRSVS,4)
CALL BUDGET (PRSVS(:,:,:,JL),JL+12,'FRC_BU_RSV')
END DO
END IF
!
!----------------------------------------------------------------------------
!
! deallocate work arrays
!
DEALLOCATE(ZUF)
DEALLOCATE(ZVF)
DEALLOCATE(ZWF)
DEALLOCATE(ZTHF)
DEALLOCATE(ZRVF)
DEALLOCATE(ZGXTHF)
DEALLOCATE(ZGYTHF)
DEALLOCATE(ZTENDTHF)
DEALLOCATE(ZTENDRVF)
DEALLOCATE(ZDZZ)
DEALLOCATE(ZRWCF)
DEALLOCATE(ZDUF)
DEALLOCATE(ZDVF)
!
!----------------------------------------------------------------------------
!
CONTAINS
SUBROUTINE DEFINE_RELAX_FORCING( OMASK_RELAX_FRC,PTHT,PDZZ, &
KPHEXT,KPVEXT )
!
LOGICAL, DIMENSION(:,:,:), INTENT(OUT) :: OMASK_RELAX_FRC
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT
REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ
INTEGER, INTENT(IN) :: KPHEXT, KPVEXT
!
INTEGER :: JI, JJ, JK, IIBEG, IIEND, IJBEG, IJEND, IKBEG, IKEND
INTEGER :: IKGRAD_TH_MAX
REAL :: ZGRAD_TH, ZGRAD_TH_MAX
!
IIBEG = KPHEXT + 1
IIEND = SIZE(PTHT,1) - KPHEXT
IJBEG = KPHEXT + 1
IJEND = SIZE(PTHT,2) - KPHEXT
IKBEG = KPHEXT + 1
IKEND = SIZE(PTHT,3) - KPVEXT
!
OMASK_RELAX_FRC(:,:,:) = .TRUE.
!
DO JI = IIBEG , IIEND
DO JJ = IJBEG , IJEND
ZGRAD_TH_MAX = -1.E30
DO JK = IKBEG+1 , IKEND-1
ZGRAD_TH = (PTHT(JI,JJ,JK+1)-PTHT(JI,JJ,JK))/PDZZ(JI,JJ,JK+1)
IF( ZGRAD_TH > ZGRAD_TH_MAX ) THEN
IKGRAD_TH_MAX = JK
ZGRAD_TH_MAX = ZGRAD_TH
END IF
END DO
IKGRAD_TH_MAX = MIN( SIZE(PTHT,3),IKGRAD_TH_MAX+2 )
OMASK_RELAX_FRC(:,:,1:IKGRAD_TH_MAX) = .FALSE.
END DO
END DO
!
IF (NVERB >= 10) THEN
CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT0,IRESP)
WRITE(ILUOUT0,*) 'DEFINE_RELAX_FORCING: IKGRAD_TH_MAX = ',IKGRAD_TH_MAX
END IF
!
END SUBROUTINE DEFINE_RELAX_FORCING
!
!----------------------------------------------------------------------------
!
END SUBROUTINE FORCING