source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/yamada4.F90 @ 3820

Last change on this file since 3820 was 3818, checked in by millour, 10 years ago

Some partial cleanup on uses of "dimensions.h" in physics.
At this point 3D gcm compiles and bench seems to run fine :-)
EM

File size: 18.1 KB
Line 
1
2! $Header$
3
4SUBROUTINE yamada4(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
5    cd, q2, km, kn, kq, ustar, iflag_pbl)
6  USE dimphy
7  USE inifis_mod, ONLY: prt_level
8  IMPLICIT NONE
9
10  ! dt : pas de temps
11  ! g  : g
12  ! zlev : altitude a chaque niveau (interface inferieure de la couche
13  ! de meme indice)
14  ! zlay : altitude au centre de chaque couche
15  ! u,v : vitesse au centre de chaque couche
16  ! (en entree : la valeur au debut du pas de temps)
17  ! teta : temperature potentielle au centre de chaque couche
18  ! (en entree : la valeur au debut du pas de temps)
19  ! cd : cdrag
20  ! (en entree : la valeur au debut du pas de temps)
21  ! q2 : $q^2$ au bas de chaque couche
22  ! (en entree : la valeur au debut du pas de temps)
23  ! (en sortie : la valeur a la fin du pas de temps)
24  ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
25  ! couche)
26  ! (en sortie : la valeur a la fin du pas de temps)
27  ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
28  ! (en sortie : la valeur a la fin du pas de temps)
29
30  ! iflag_pbl doit valoir entre 6 et 9
31  ! l=6, on prend  systematiquement une longueur d'equilibre
32  ! iflag_pbl=6 : MY 2.0
33  ! iflag_pbl=7 : MY 2.0.Fournier
34  ! iflag_pbl=8/9 : MY 2.5
35  ! iflag_pbl=8 with special obsolete treatments for convergence
36  ! with Cmpi5 NPv3.1 simulations
37  ! iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact
38  ! iflag_pbl=12 = 11 with vertical diffusion off q2
39
40  ! 2013/04/01 (FH hourdin@lmd.jussieu.fr)
41  ! Correction for very stable PBLs (iflag_pbl=10 and 11)
42  ! iflag_pbl=8 converges numerically with NPv3.1
43  ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l
44  ! -> the model can run with longer time-steps.
45  ! .......................................................................
46
47  REAL dt, g, rconst
48  REAL plev(klon, klev+1), temp(klon, klev)
49  REAL ustar(klon)
50  REAL kmin, qmin, pblhmin(klon), coriol(klon)
51  REAL zlev(klon, klev+1)
52  REAL zlay(klon, klev)
53  REAL u(klon, klev)
54  REAL v(klon, klev)
55  REAL teta(klon, klev)
56  REAL cd(klon)
57  REAL q2(klon, klev+1), qpre
58  REAL unsdz(klon, klev)
59  REAL unsdzdec(klon, klev+1)
60
61  REAL km(klon, klev+1)
62  REAL kmpre(klon, klev+1), tmp2
63  REAL mpre(klon, klev+1)
64  REAL kn(klon, klev+1)
65  REAL kq(klon, klev+1)
66  REAL ff(klon, klev+1), delta(klon, klev+1)
67  REAL aa(klon, klev+1), aa0, aa1
68  INTEGER iflag_pbl, ngrid
69  INTEGER nlay, nlev
70
71  LOGICAL first
72  INTEGER ipas
73  SAVE first, ipas
74  ! FH/IM     data first,ipas/.true.,0/
75  DATA first, ipas/.FALSE., 0/
76  !$OMP THREADPRIVATE( first,ipas)
77
78  INTEGER ig, k
79
80
81  REAL ri, zrif, zalpha, zsm, zsn
82  REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
83
84  REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
85  REAL dtetadz(klon, klev+1)
86  REAL m2cstat, mcstat, kmcstat
87  REAL l(klon, klev+1)
88  REAL, ALLOCATABLE, SAVE :: l0(:)
89  !$OMP THREADPRIVATE(l0)
90  REAL sq(klon), sqz(klon), zz(klon, klev+1)
91  INTEGER iter
92
93  REAL ric, rifc, b1, kap
94  SAVE ric, rifc, b1, kap
95  DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/
96  !$OMP THREADPRIVATE(ric,rifc,b1,kap)
97  REAL frif, falpha, fsm
98  REAL fl, zzz, zl0, zq2, zn2
99
100  REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &
101    lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev)
102  LOGICAL, SAVE :: firstcall = .TRUE.
103  !$OMP THREADPRIVATE(firstcall)
104  frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
105  falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
106  fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
107  fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, &
108    k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( &
109    max(n2(ig,k),1.E-10))), 1.)
110
111
112  nlay = klev
113  nlev = klev + 1
114
115  IF (firstcall) THEN
116    ALLOCATE (l0(klon))
117    firstcall = .FALSE.
118  END IF
119
120
121  IF (.NOT. (iflag_pbl>=6 .AND. iflag_pbl<=12)) THEN
122    STOP 'probleme de coherence dans appel a MY'
123  END IF
124
125  ipas = ipas + 1
126
127
128  ! .......................................................................
129  ! les increments verticaux
130  ! .......................................................................
131
132  ! !!!!! allerte !!!!!c
133  ! !!!!! zlev n'est pas declare a nlev !!!!!c
134  ! !!!!! ---->
135  DO ig = 1, ngrid
136    zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1))
137  END DO
138  ! !!!!! <----
139  ! !!!!! allerte !!!!!c
140
141  DO k = 1, nlay
142    DO ig = 1, ngrid
143      unsdz(ig, k) = 1.E+0/(zlev(ig,k+1)-zlev(ig,k))
144    END DO
145  END DO
146  DO ig = 1, ngrid
147    unsdzdec(ig, 1) = 1.E+0/(zlay(ig,1)-zlev(ig,1))
148  END DO
149  DO k = 2, nlay
150    DO ig = 1, ngrid
151      unsdzdec(ig, k) = 1.E+0/(zlay(ig,k)-zlay(ig,k-1))
152    END DO
153  END DO
154  DO ig = 1, ngrid
155    unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
156  END DO
157
158  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
159  ! Computing M^2, N^2, Richardson numbers, stability functions
160  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
161    ! initialize arrays:
162  m2(:, :) = 0.0
163  sm(:, :) = 0.0
164  rif(:, :) = 0.0
165
166  DO k = 2, klev
167    DO ig = 1, ngrid
168      dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
169      m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &
170        k-1))**2)/(dz(ig,k)*dz(ig,k))
171      dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k)
172      n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k))
173      ! n2(ig,k)=0.
174      ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
175      IF (ri<ric) THEN
176        rif(ig, k) = frif(ri)
177      ELSE
178        rif(ig, k) = rifc
179      END IF
180      IF (rif(ig,k)<0.16) THEN
181        alpha(ig, k) = falpha(rif(ig,k))
182        sm(ig, k) = fsm(rif(ig,k))
183      ELSE
184        alpha(ig, k) = 1.12
185        sm(ig, k) = 0.085
186      END IF
187      zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
188    END DO
189  END DO
190
191
192  ! ====================================================================
193  ! Computing the mixing length
194  ! ====================================================================
195
196  ! Mise a jour de l0
197  IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
198
199    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
200    ! Iterative computation of l0
201    ! This version is kept for iflag_pbl only for convergence
202    ! with NPv3.1 Cmip5 simulations
203    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
204
205    DO ig = 1, ngrid
206      sq(ig) = 1.E-10
207      sqz(ig) = 1.E-10
208    END DO
209    DO k = 2, klev - 1
210      DO ig = 1, ngrid
211        zq = sqrt(q2(ig,k))
212        sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
213        sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
214      END DO
215    END DO
216    DO ig = 1, ngrid
217      l0(ig) = 0.2*sqz(ig)/sq(ig)
218    END DO
219    DO k = 2, klev
220      DO ig = 1, ngrid
221        l(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
222      END DO
223    END DO
224    ! print*,'L0 cas 8 ou 10 ',l0
225
226  ELSE
227
228    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
229    ! In all other case, the assymptotic mixing length l0 is imposed (100m)
230    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
231
232    l0(:) = 150.
233    DO k = 2, klev
234      DO ig = 1, ngrid
235        l(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
236      END DO
237    END DO
238    ! print*,'L0 cas autres ',l0
239
240  END IF
241
242
243  ! ====================================================================
244  ! Yamada 2.0
245  ! ====================================================================
246  IF (iflag_pbl==6) THEN
247
248    DO k = 2, klev
249      q2(:, k) = l(:, k)**2*zz(:, k)
250    END DO
251
252
253  ELSE IF (iflag_pbl==7) THEN
254    ! ====================================================================
255    ! Yamada 2.Fournier
256    ! ====================================================================
257
258    ! Calcul de l,  km, au pas precedent
259    DO k = 2, klev
260      DO ig = 1, ngrid
261        ! print*,'SMML=',sm(ig,k),l(ig,k)
262        delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
263        kmpre(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
264        mpre(ig, k) = sqrt(m2(ig,k))
265        ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
266      END DO
267    END DO
268
269    DO k = 2, klev - 1
270      DO ig = 1, ngrid
271        m2cstat = max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1, 1.E-12)
272        mcstat = sqrt(m2cstat)
273
274        ! print*,'M2 L=',k,mpre(ig,k),mcstat
275
276        ! -----{puis on ecrit la valeur de q qui annule l'equation de m
277        ! supposee en q3}
278
279        IF (k==2) THEN
280          kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ &
281            unsdz(ig,k-1)*cd(ig)*(sqrt(u(ig,3)**2+v(ig,3)**2)-mcstat/unsdzdec &
282            (ig,k)-mpre(ig,k+1)/unsdzdec(ig,k+1))**2)/(unsdz(ig,k)+unsdz(ig,k &
283            -1))
284        ELSE
285          kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ &
286            unsdz(ig,k-1)*kmpre(ig,k-1)*mpre(ig,k-1))/ &
287            (unsdz(ig,k)+unsdz(ig,k-1))
288        END IF
289        ! print*,'T2 L=',k,tmp2
290        tmp2 = kmcstat/(sm(ig,k)/q2(ig,k))/l(ig, k)
291        q2(ig, k) = max(tmp2, 1.E-12)**(2./3.)
292        ! print*,'Q2 L=',k,q2(ig,k)
293
294      END DO
295    END DO
296
297  ELSE IF (iflag_pbl==8 .OR. iflag_pbl==9) THEN
298    ! ====================================================================
299    ! Yamada 2.5 a la Didi
300    ! ====================================================================
301
302
303    ! Calcul de l,  km, au pas precedent
304    DO k = 2, klev
305      DO ig = 1, ngrid
306        ! print*,'SMML=',sm(ig,k),l(ig,k)
307        delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
308        IF (delta(ig,k)<1.E-20) THEN
309          ! print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
310          delta(ig, k) = 1.E-20
311        END IF
312        km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
313        aa0 = (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
314        aa1 = (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
315        ! abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
316        aa(ig, k) = aa1*dt/(delta(ig,k)*l(ig,k))
317        ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
318        qpre = sqrt(q2(ig,k))
319        ! if (iflag_pbl.eq.8 ) then
320        IF (aa(ig,k)>0.) THEN
321          q2(ig, k) = (qpre+aa(ig,k)*qpre*qpre)**2
322        ELSE
323          q2(ig, k) = (qpre/(1.-aa(ig,k)*qpre))**2
324        END IF
325        ! else ! iflag_pbl=9
326        ! if (aa(ig,k)*qpre.gt.0.9) then
327        ! q2(ig,k)=(qpre*10.)**2
328        ! else
329        ! q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
330        ! endif
331        ! endif
332        q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
333        ! print*,'Q2 L=',k,q2(ig,k),qpre*qpre
334      END DO
335    END DO
336
337  ELSE IF (iflag_pbl>=10) THEN
338
339    ! print*,'Schema mixte D'
340    ! print*,'Longueur ',l(:,:)
341    DO k = 2, klev - 1
342      l(:, k) = max(l(:,k), 1.)
343      km(:, k) = l(:, k)*sqrt(q2(:,k))*sm(:, k)
344      q2(:, k) = q2(:, k) + dt*km(:, k)*m2(:, k)*(1.-rif(:,k))
345      q2(:, k) = min(max(q2(:,k),1.E-10), 1.E4)
346      q2(:, k) = 1./(1./sqrt(q2(:,k))+dt/(2*l(:,k)*b1))
347      q2(:, k) = q2(:, k)*q2(:, k)
348    END DO
349
350
351  ELSE
352    STOP 'Cas nom prevu dans yamada4'
353
354  END IF ! Fin du cas 8
355
356  ! print*,'OK8'
357
358  ! ====================================================================
359  ! Calcul des coefficients de m�ange
360  ! ====================================================================
361  DO k = 2, klev
362    ! print*,'k=',k
363    DO ig = 1, ngrid
364      ! abde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
365      zq = sqrt(q2(ig,k))
366      km(ig, k) = l(ig, k)*zq*sm(ig, k)
367      kn(ig, k) = km(ig, k)*alpha(ig, k)
368      kq(ig, k) = l(ig, k)*zq*0.2
369      ! print*,'KML=',km(ig,k),kn(ig,k)
370    END DO
371  END DO
372    ! initialize near-surface and top-layer mixing coefficients
373  kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface
374  kq(1:ngrid, klev+1) = 0 ! zero at the top
375
376  ! Transport diffusif vertical de la TKE.
377  IF (iflag_pbl>=12) THEN
378    ! print*,'YAMADA VDIF'
379    q2(:, 1) = q2(:, 2)
380    CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2)
381  END IF
382
383  ! Traitement des cas noctrunes avec l'introduction d'une longueur
384  ! minilale.
385
386  ! ====================================================================
387  ! Traitement particulier pour les cas tres stables.
388  ! D'apres Holtslag Boville.
389
390  IF (prt_level>1) THEN
391    PRINT *, 'YAMADA4 0'
392  END IF !(prt_level>1) THEN
393  DO ig = 1, ngrid
394    coriol(ig) = 1.E-4
395    pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5)
396  END DO
397
398  ! print*,'pblhmin ',pblhmin
399  ! Test a remettre 21 11 02
400  ! test abd 13 05 02      if(0.eq.1) then
401  IF (1==1) THEN
402    IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
403
404      DO k = 2, klev
405        DO ig = 1, ngrid
406          IF (teta(ig,2)>teta(ig,1)) THEN
407            qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
408            kmin = kap*zlev(ig, k)*qmin
409          ELSE
410            kmin = -1. ! kmin n'est utilise que pour les SL stables.
411          END IF
412          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
413            ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
414            ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
415            kn(ig, k) = kmin
416            km(ig, k) = kmin
417            kq(ig, k) = kmin
418            ! la longueur de melange est suposee etre l= kap z
419            ! K=l q Sm d'ou q2=(K/l Sm)**2
420            q2(ig, k) = (qmin/sm(ig,k))**2
421          END IF
422        END DO
423      END DO
424
425    ELSE
426
427      DO k = 2, klev
428        DO ig = 1, ngrid
429          IF (teta(ig,2)>teta(ig,1)) THEN
430            qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
431            kmin = kap*zlev(ig, k)*qmin
432          ELSE
433            kmin = -1. ! kmin n'est utilise que pour les SL stables.
434          END IF
435          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
436            ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
437            ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
438            kn(ig, k) = kmin
439            km(ig, k) = kmin
440            kq(ig, k) = kmin
441            ! la longueur de melange est suposee etre l= kap z
442            ! K=l q Sm d'ou q2=(K/l Sm)**2
443            sm(ig, k) = 1.
444            alpha(ig, k) = 1.
445            q2(ig, k) = min((qmin/sm(ig,k))**2, 10.)
446            zq = sqrt(q2(ig,k))
447            km(ig, k) = l(ig, k)*zq*sm(ig, k)
448            kn(ig, k) = km(ig, k)*alpha(ig, k)
449            kq(ig, k) = l(ig, k)*zq*0.2
450          END IF
451        END DO
452      END DO
453    END IF
454
455  END IF
456
457  IF (prt_level>1) THEN
458    PRINT *, 'YAMADA4 1'
459  END IF !(prt_level>1) THEN
460  ! Diagnostique pour stokage
461
462  IF (1==0) THEN
463    rino = rif
464    smyam(1:ngrid, 1) = 0.
465    styam(1:ngrid, 1) = 0.
466    lyam(1:ngrid, 1) = 0.
467    knyam(1:ngrid, 1) = 0.
468    w2yam(1:ngrid, 1) = 0.
469    t2yam(1:ngrid, 1) = 0.
470
471    smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)
472    styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev)
473    lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev)
474    knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)
475
476    ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
477
478    w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &
479      lyam(1:ngrid, 2:klev)*5.17*kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ &
480      sqrt(q2(1:ngrid,2:klev))
481
482    t2yam(1:ngrid, 2:klev) = 9.1*kn(1:ngrid, 2:klev)* &
483      dtetadz(1:ngrid, 2:klev)**2/sqrt(q2(1:ngrid,2:klev))* &
484      lyam(1:ngrid, 2:klev)
485  END IF
486
487  ! print*,'OKFIN'
488  first = .FALSE.
489  RETURN
490END SUBROUTINE yamada4
491SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
492  USE dimphy
493  IMPLICIT NONE
494
495  ! dt : pas de temps
496
497  REAL plev(klon, klev+1)
498  REAL temp(klon, klev)
499  REAL timestep
500  REAL gravity, rconst
501  REAL kstar(klon, klev+1), zz
502  REAL kmy(klon, klev+1)
503  REAL q2(klon, klev+1)
504  REAL deltap(klon, klev+1)
505  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
506  INTEGER ngrid
507
508  INTEGER i, k
509
510  ! print*,'RD=',rconst
511  DO k = 1, klev
512    DO i = 1, ngrid
513      ! test
514      ! print*,'i,k',i,k
515      ! print*,'temp(i,k)=',temp(i,k)
516      ! print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
517      zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
518      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
519        (plev(i,k)-plev(i,k+1))*timestep
520    END DO
521  END DO
522
523  DO k = 2, klev
524    DO i = 1, ngrid
525      deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
526    END DO
527  END DO
528  DO i = 1, ngrid
529    deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
530    deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
531    denom(i, klev+1) = deltap(i, klev+1) + kstar(i, klev)
532    alpha(i, klev+1) = deltap(i, klev+1)*q2(i, klev+1)/denom(i, klev+1)
533    beta(i, klev+1) = kstar(i, klev)/denom(i, klev+1)
534  END DO
535
536  DO k = klev, 2, -1
537    DO i = 1, ngrid
538      denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + &
539        kstar(i, k-1)
540      ! correction d'un bug 10 01 2001
541      alpha(i, k) = (q2(i,k)*deltap(i,k)+kstar(i,k)*alpha(i,k+1))/denom(i, k)
542      beta(i, k) = kstar(i, k-1)/denom(i, k)
543    END DO
544  END DO
545
546  ! Si on recalcule q2(1)
547  IF (1==0) THEN
548    DO i = 1, ngrid
549      denom(i, 1) = deltap(i, 1) + (1-beta(i,2))*kstar(i, 1)
550      q2(i, 1) = (q2(i,1)*deltap(i,1)+kstar(i,1)*alpha(i,2))/denom(i, 1)
551    END DO
552  END IF
553  ! sinon, on peut sauter cette boucle...
554
555  DO k = 2, klev + 1
556    DO i = 1, ngrid
557      q2(i, k) = alpha(i, k) + beta(i, k)*q2(i, k-1)
558    END DO
559  END DO
560
561  RETURN
562END SUBROUTINE vdif_q2
563SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
564  USE dimphy
565  IMPLICIT NONE
566
567  ! dt : pas de temps
568
569  REAL plev(klon, klev+1)
570  REAL temp(klon, klev)
571  REAL timestep
572  REAL gravity, rconst
573  REAL kstar(klon, klev+1), zz
574  REAL kmy(klon, klev+1)
575  REAL q2(klon, klev+1)
576  REAL deltap(klon, klev+1)
577  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
578  INTEGER ngrid
579
580  INTEGER i, k
581
582  DO k = 1, klev
583    DO i = 1, ngrid
584      zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
585      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
586        (plev(i,k)-plev(i,k+1))*timestep
587    END DO
588  END DO
589
590  DO k = 2, klev
591    DO i = 1, ngrid
592      deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
593    END DO
594  END DO
595  DO i = 1, ngrid
596    deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
597    deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
598  END DO
599
600  DO k = klev, 2, -1
601    DO i = 1, ngrid
602      q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, &
603        k))-kstar(i,k-1)*(q2(i,k)-q2(i,k-1)))/deltap(i, k)
604    END DO
605  END DO
606
607  DO i = 1, ngrid
608    q2(i, 1) = q2(i, 1) + (kstar(i,1)*(q2(i,2)-q2(i,1)))/deltap(i, 1)
609    q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, &
610      klev)))/deltap(i, klev+1)
611  END DO
612
613  RETURN
614END SUBROUTINE vdif_q2e
Note: See TracBrowser for help on using the repository browser.