source: LMDZ5/branches/testing/libf/phylmd/yamada4.F90 @ 2408

Last change on this file since 2408 was 2408, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2298:2396 into testing branch
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
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 print_control_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      DO ig = 1, ngrid
343        l(ig, k) = max(l(ig,k), 1.)
344        km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
345        q2(ig, k) = q2(ig, k) + dt*km(ig, k)*m2(ig, k)*(1.-rif(ig,k))
346        q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
347        q2(ig, k) = 1./(1./sqrt(q2(ig,k))+dt/(2*l(ig,k)*b1))
348        q2(ig, k) = q2(ig, k)*q2(ig, k)
349      END DO
350    END DO
351
352
353  ELSE
354    STOP 'Cas nom prevu dans yamada4'
355
356  END IF ! Fin du cas 8
357
358  ! print*,'OK8'
359
360  ! ====================================================================
361  ! Calcul des coefficients de m�ange
362  ! ====================================================================
363  DO k = 2, klev
364    ! print*,'k=',k
365    DO ig = 1, ngrid
366      ! abde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
367      zq = sqrt(q2(ig,k))
368      km(ig, k) = l(ig, k)*zq*sm(ig, k)
369      kn(ig, k) = km(ig, k)*alpha(ig, k)
370      kq(ig, k) = l(ig, k)*zq*0.2
371      ! print*,'KML=',km(ig,k),kn(ig,k)
372    END DO
373  END DO
374    ! initialize near-surface and top-layer mixing coefficients
375  kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface
376  kq(1:ngrid, klev+1) = 0 ! zero at the top
377
378  ! Transport diffusif vertical de la TKE.
379  IF (iflag_pbl>=12) THEN
380    ! print*,'YAMADA VDIF'
381    q2(:, 1) = q2(:, 2)
382    CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2)
383  END IF
384
385  ! Traitement des cas noctrunes avec l'introduction d'une longueur
386  ! minilale.
387
388  ! ====================================================================
389  ! Traitement particulier pour les cas tres stables.
390  ! D'apres Holtslag Boville.
391
392  IF (prt_level>1) THEN
393    PRINT *, 'YAMADA4 0'
394  END IF !(prt_level>1) THEN
395  DO ig = 1, ngrid
396    coriol(ig) = 1.E-4
397    pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5)
398  END DO
399
400  ! print*,'pblhmin ',pblhmin
401  ! Test a remettre 21 11 02
402  ! test abd 13 05 02      if(0.eq.1) then
403  IF (1==1) THEN
404    IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
405
406      DO k = 2, klev
407        DO ig = 1, ngrid
408          IF (teta(ig,2)>teta(ig,1)) THEN
409            qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
410            kmin = kap*zlev(ig, k)*qmin
411          ELSE
412            kmin = -1. ! kmin n'est utilise que pour les SL stables.
413          END IF
414          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
415            ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
416            ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
417            kn(ig, k) = kmin
418            km(ig, k) = kmin
419            kq(ig, k) = kmin
420            ! la longueur de melange est suposee etre l= kap z
421            ! K=l q Sm d'ou q2=(K/l Sm)**2
422            q2(ig, k) = (qmin/sm(ig,k))**2
423          END IF
424        END DO
425      END DO
426
427    ELSE
428
429      DO k = 2, klev
430        DO ig = 1, ngrid
431          IF (teta(ig,2)>teta(ig,1)) THEN
432            qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
433            kmin = kap*zlev(ig, k)*qmin
434          ELSE
435            kmin = -1. ! kmin n'est utilise que pour les SL stables.
436          END IF
437          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
438            ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
439            ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
440            kn(ig, k) = kmin
441            km(ig, k) = kmin
442            kq(ig, k) = kmin
443            ! la longueur de melange est suposee etre l= kap z
444            ! K=l q Sm d'ou q2=(K/l Sm)**2
445            sm(ig, k) = 1.
446            alpha(ig, k) = 1.
447            q2(ig, k) = min((qmin/sm(ig,k))**2, 10.)
448            zq = sqrt(q2(ig,k))
449            km(ig, k) = l(ig, k)*zq*sm(ig, k)
450            kn(ig, k) = km(ig, k)*alpha(ig, k)
451            kq(ig, k) = l(ig, k)*zq*0.2
452          END IF
453        END DO
454      END DO
455    END IF
456
457  END IF
458
459  IF (prt_level>1) THEN
460    PRINT *, 'YAMADA4 1'
461  END IF !(prt_level>1) THEN
462  ! Diagnostique pour stokage
463
464  IF (1==0) THEN
465    rino = rif
466    smyam(1:ngrid, 1) = 0.
467    styam(1:ngrid, 1) = 0.
468    lyam(1:ngrid, 1) = 0.
469    knyam(1:ngrid, 1) = 0.
470    w2yam(1:ngrid, 1) = 0.
471    t2yam(1:ngrid, 1) = 0.
472
473    smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)
474    styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev)
475    lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev)
476    knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)
477
478    ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
479
480    w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &
481      lyam(1:ngrid, 2:klev)*5.17*kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ &
482      sqrt(q2(1:ngrid,2:klev))
483
484    t2yam(1:ngrid, 2:klev) = 9.1*kn(1:ngrid, 2:klev)* &
485      dtetadz(1:ngrid, 2:klev)**2/sqrt(q2(1:ngrid,2:klev))* &
486      lyam(1:ngrid, 2:klev)
487  END IF
488
489  ! print*,'OKFIN'
490  first = .FALSE.
491  RETURN
492END SUBROUTINE yamada4
493SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
494  USE dimphy
495  IMPLICIT NONE
496
497  ! dt : pas de temps
498
499  REAL plev(klon, klev+1)
500  REAL temp(klon, klev)
501  REAL timestep
502  REAL gravity, rconst
503  REAL kstar(klon, klev+1), zz
504  REAL kmy(klon, klev+1)
505  REAL q2(klon, klev+1)
506  REAL deltap(klon, klev+1)
507  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
508  INTEGER ngrid
509
510  INTEGER i, k
511
512  ! print*,'RD=',rconst
513  DO k = 1, klev
514    DO i = 1, ngrid
515      ! test
516      ! print*,'i,k',i,k
517      ! print*,'temp(i,k)=',temp(i,k)
518      ! print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
519      zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
520      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
521        (plev(i,k)-plev(i,k+1))*timestep
522    END DO
523  END DO
524
525  DO k = 2, klev
526    DO i = 1, ngrid
527      deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
528    END DO
529  END DO
530  DO i = 1, ngrid
531    deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
532    deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
533    denom(i, klev+1) = deltap(i, klev+1) + kstar(i, klev)
534    alpha(i, klev+1) = deltap(i, klev+1)*q2(i, klev+1)/denom(i, klev+1)
535    beta(i, klev+1) = kstar(i, klev)/denom(i, klev+1)
536  END DO
537
538  DO k = klev, 2, -1
539    DO i = 1, ngrid
540      denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + &
541        kstar(i, k-1)
542      ! correction d'un bug 10 01 2001
543      alpha(i, k) = (q2(i,k)*deltap(i,k)+kstar(i,k)*alpha(i,k+1))/denom(i, k)
544      beta(i, k) = kstar(i, k-1)/denom(i, k)
545    END DO
546  END DO
547
548  ! Si on recalcule q2(1)
549  IF (1==0) THEN
550    DO i = 1, ngrid
551      denom(i, 1) = deltap(i, 1) + (1-beta(i,2))*kstar(i, 1)
552      q2(i, 1) = (q2(i,1)*deltap(i,1)+kstar(i,1)*alpha(i,2))/denom(i, 1)
553    END DO
554  END IF
555  ! sinon, on peut sauter cette boucle...
556
557  DO k = 2, klev + 1
558    DO i = 1, ngrid
559      q2(i, k) = alpha(i, k) + beta(i, k)*q2(i, k-1)
560    END DO
561  END DO
562
563  RETURN
564END SUBROUTINE vdif_q2
565SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
566  USE dimphy
567  IMPLICIT NONE
568
569  ! dt : pas de temps
570
571  REAL plev(klon, klev+1)
572  REAL temp(klon, klev)
573  REAL timestep
574  REAL gravity, rconst
575  REAL kstar(klon, klev+1), zz
576  REAL kmy(klon, klev+1)
577  REAL q2(klon, klev+1)
578  REAL deltap(klon, klev+1)
579  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
580  INTEGER ngrid
581
582  INTEGER i, k
583
584  DO k = 1, klev
585    DO i = 1, ngrid
586      zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
587      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
588        (plev(i,k)-plev(i,k+1))*timestep
589    END DO
590  END DO
591
592  DO k = 2, klev
593    DO i = 1, ngrid
594      deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
595    END DO
596  END DO
597  DO i = 1, ngrid
598    deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
599    deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
600  END DO
601
602  DO k = klev, 2, -1
603    DO i = 1, ngrid
604      q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, &
605        k))-kstar(i,k-1)*(q2(i,k)-q2(i,k-1)))/deltap(i, k)
606    END DO
607  END DO
608
609  DO i = 1, ngrid
610    q2(i, 1) = q2(i, 1) + (kstar(i,1)*(q2(i,2)-q2(i,1)))/deltap(i, 1)
611    q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, &
612      klev)))/deltap(i, klev+1)
613  END DO
614
615  RETURN
616END SUBROUTINE vdif_q2e
Note: See TracBrowser for help on using the repository browser.