source: LMDZ5/trunk/libf/phylmd/yamada4.F90 @ 2345

Last change on this file since 2345 was 2339, checked in by musat, 9 years ago

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