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

Last change on this file since 2311 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

  • moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
  • created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
  • created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
  • reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

EM
  • 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      l(:, k) = max(l(:,k), 1.)
347      km(:, k) = l(:, k)*sqrt(q2(:,k))*sm(:, k)
348      q2(:, k) = q2(:, k) + dt*km(:, k)*m2(:, k)*(1.-rif(:,k))
349      q2(:, k) = min(max(q2(:,k),1.E-10), 1.E4)
350      q2(:, k) = 1./(1./sqrt(q2(:,k))+dt/(2*l(:,k)*b1))
351      q2(:, k) = q2(:, k)*q2(:, k)
352    END DO
353
354
355  ELSE
356    STOP 'Cas nom prevu dans yamada4'
357
358  END IF ! Fin du cas 8
359
360  ! print*,'OK8'
361
362  ! ====================================================================
363  ! Calcul des coefficients de m�ange
364  ! ====================================================================
365  DO k = 2, klev
366    ! print*,'k=',k
367    DO ig = 1, ngrid
368      ! abde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
369      zq = sqrt(q2(ig,k))
370      km(ig, k) = l(ig, k)*zq*sm(ig, k)
371      kn(ig, k) = km(ig, k)*alpha(ig, k)
372      kq(ig, k) = l(ig, k)*zq*0.2
373      ! print*,'KML=',km(ig,k),kn(ig,k)
374    END DO
375  END DO
376    ! initialize near-surface and top-layer mixing coefficients
377  kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface
378  kq(1:ngrid, klev+1) = 0 ! zero at the top
379
380  ! Transport diffusif vertical de la TKE.
381  IF (iflag_pbl>=12) THEN
382    ! print*,'YAMADA VDIF'
383    q2(:, 1) = q2(:, 2)
384    CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2)
385  END IF
386
387  ! Traitement des cas noctrunes avec l'introduction d'une longueur
388  ! minilale.
389
390  ! ====================================================================
391  ! Traitement particulier pour les cas tres stables.
392  ! D'apres Holtslag Boville.
393
394  IF (prt_level>1) THEN
395    PRINT *, 'YAMADA4 0'
396  END IF !(prt_level>1) THEN
397  DO ig = 1, ngrid
398    coriol(ig) = 1.E-4
399    pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5)
400  END DO
401
402  ! print*,'pblhmin ',pblhmin
403  ! Test a remettre 21 11 02
404  ! test abd 13 05 02      if(0.eq.1) then
405  IF (1==1) THEN
406    IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
407
408      DO k = 2, klev
409        DO ig = 1, ngrid
410          IF (teta(ig,2)>teta(ig,1)) THEN
411            qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
412            kmin = kap*zlev(ig, k)*qmin
413          ELSE
414            kmin = -1. ! kmin n'est utilise que pour les SL stables.
415          END IF
416          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
417            ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
418            ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
419            kn(ig, k) = kmin
420            km(ig, k) = kmin
421            kq(ig, k) = kmin
422            ! la longueur de melange est suposee etre l= kap z
423            ! K=l q Sm d'ou q2=(K/l Sm)**2
424            q2(ig, k) = (qmin/sm(ig,k))**2
425          END IF
426        END DO
427      END DO
428
429    ELSE
430
431      DO k = 2, klev
432        DO ig = 1, ngrid
433          IF (teta(ig,2)>teta(ig,1)) THEN
434            qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
435            kmin = kap*zlev(ig, k)*qmin
436          ELSE
437            kmin = -1. ! kmin n'est utilise que pour les SL stables.
438          END IF
439          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
440            ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
441            ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
442            kn(ig, k) = kmin
443            km(ig, k) = kmin
444            kq(ig, k) = kmin
445            ! la longueur de melange est suposee etre l= kap z
446            ! K=l q Sm d'ou q2=(K/l Sm)**2
447            sm(ig, k) = 1.
448            alpha(ig, k) = 1.
449            q2(ig, k) = min((qmin/sm(ig,k))**2, 10.)
450            zq = sqrt(q2(ig,k))
451            km(ig, k) = l(ig, k)*zq*sm(ig, k)
452            kn(ig, k) = km(ig, k)*alpha(ig, k)
453            kq(ig, k) = l(ig, k)*zq*0.2
454          END IF
455        END DO
456      END DO
457    END IF
458
459  END IF
460
461  IF (prt_level>1) THEN
462    PRINT *, 'YAMADA4 1'
463  END IF !(prt_level>1) THEN
464  ! Diagnostique pour stokage
465
466  IF (1==0) THEN
467    rino = rif
468    smyam(1:ngrid, 1) = 0.
469    styam(1:ngrid, 1) = 0.
470    lyam(1:ngrid, 1) = 0.
471    knyam(1:ngrid, 1) = 0.
472    w2yam(1:ngrid, 1) = 0.
473    t2yam(1:ngrid, 1) = 0.
474
475    smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)
476    styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev)
477    lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev)
478    knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)
479
480    ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
481
482    w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &
483      lyam(1:ngrid, 2:klev)*5.17*kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ &
484      sqrt(q2(1:ngrid,2:klev))
485
486    t2yam(1:ngrid, 2:klev) = 9.1*kn(1:ngrid, 2:klev)* &
487      dtetadz(1:ngrid, 2:klev)**2/sqrt(q2(1:ngrid,2:klev))* &
488      lyam(1:ngrid, 2:klev)
489  END IF
490
491  ! print*,'OKFIN'
492  first = .FALSE.
493  RETURN
494END SUBROUTINE yamada4
495SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
496  USE dimphy
497  IMPLICIT NONE
498  ! .......................................................................
499  include "dimensions.h"
500  ! ccc#include "dimphy.h"
501  ! .......................................................................
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  include "dimensions.h"
576  ! ccc#include "dimphy.h"
577  ! .......................................................................
578
579  ! dt : pas de temps
580
581  REAL plev(klon, klev+1)
582  REAL temp(klon, klev)
583  REAL timestep
584  REAL gravity, rconst
585  REAL kstar(klon, klev+1), zz
586  REAL kmy(klon, klev+1)
587  REAL q2(klon, klev+1)
588  REAL deltap(klon, klev+1)
589  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
590  INTEGER ngrid
591
592  INTEGER i, k
593
594  DO k = 1, klev
595    DO i = 1, ngrid
596      zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
597      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
598        (plev(i,k)-plev(i,k+1))*timestep
599    END DO
600  END DO
601
602  DO k = 2, klev
603    DO i = 1, ngrid
604      deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
605    END DO
606  END DO
607  DO i = 1, ngrid
608    deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
609    deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
610  END DO
611
612  DO k = klev, 2, -1
613    DO i = 1, ngrid
614      q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, &
615        k))-kstar(i,k-1)*(q2(i,k)-q2(i,k-1)))/deltap(i, k)
616    END DO
617  END DO
618
619  DO i = 1, ngrid
620    q2(i, 1) = q2(i, 1) + (kstar(i,1)*(q2(i,2)-q2(i,1)))/deltap(i, 1)
621    q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, &
622      klev)))/deltap(i, klev+1)
623  END DO
624
625  RETURN
626END SUBROUTINE vdif_q2e
Note: See TracBrowser for help on using the repository browser.