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

Last change on this file since 2237 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.
  • 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  IMPLICIT NONE
8  include "iniprint.h"
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.