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

Last change on this file since 2471 was 2471, checked in by Laurent Fairhead, 8 years ago

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