source: LMDZ5/trunk/libf/phylmd/coef_diff_turb_mod.F90 @ 1744

Last change on this file since 1744 was 1738, checked in by idelkadi, 12 years ago

Modifications for numerical stability of the boundary layer parameterizations.
Concerning yamada4.F :

  • option with asymptotic mixing length l0 imposed (iflag_pbl=8 and 9)
  • option with a new temporal scheme (iflag_pbl=10 and 11)
  • Correction for very stable PBLs (iflag_pbl=10 and 11) iflag_pbl=8 converges numerically with NPv3.1 iflag_pbl=11 -> now starts with NP from start files created by ce0l

-> the model can run with longer time-steps.

Concerning thermals :

Introduction of an implicit computation of vertical advection in
the environment of thermal plumes in thermcell_dq
impl = 0 : explicit, 1 : implicit, -1 : old version
controled by iflag_thermals =

15, 16 run with impl=-1 : numerical convergence with NPv3
17, 18 run with impl=1 : more stable

15 and 17 correspond to the activation of the stratocumulus "bidouille"

Modified routines (phylmd/):
calltherm.F90 : for managing the various options of thermcell_dq
coef_diff_turb_mod.F90 : yamada4 called for iflag_pbl<= 18 instead of 11
physiq.F : desactivation of the vertical diffusion of TKE for iflag_pbl=10
thermcellV0_main.F90 : calling thermcell_dq with implicit scheme
thermcell_dq.F90 : thermcell_dq with optional explicit/implicit scheme
thermcell_main.F90 : various option for vertical transport by thermals
yamada4.F : Yamada scheme with new options

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 18.8 KB
Line 
1!
2MODULE coef_diff_turb_mod
3!
4! This module contains some procedures for calculation of the coefficients of the
5! turbulent diffusion in the atmosphere and coefficients for turbulent diffusion
6! at surface(cdrag)
7!
8  IMPLICIT NONE
9 
10CONTAINS
11!
12!****************************************************************************************
13!
14  SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, &
15       ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
16       ycoefm, ycoefh ,yq2)
17 
18    USE dimphy
19!
20! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
21! atmosphere
22! NB! No values are calculated between surface and the first model layer.
23!     ycoefm(:,1) and ycoefh(:,1) are not valid !!!
24!
25!
26! Input arguments
27!****************************************************************************************
28    REAL, INTENT(IN)                           :: dtime
29    INTEGER, INTENT(IN)                        :: nsrf, knon
30    INTEGER, DIMENSION(klon), INTENT(IN)       :: ni
31    REAL, DIMENSION(klon,klev+1), INTENT(IN)   :: ypaprs
32    REAL, DIMENSION(klon,klev), INTENT(IN)     :: ypplay
33    REAL, DIMENSION(klon,klev), INTENT(IN)     :: yu, yv
34    REAL, DIMENSION(klon,klev), INTENT(IN)     :: yq, yt
35    REAL, DIMENSION(klon), INTENT(IN)          :: yts, yrugos, yqsurf
36    REAL, DIMENSION(klon), INTENT(IN)          :: ycdragm
37
38! InOutput arguments
39!****************************************************************************************
40    REAL, DIMENSION(klon,klev+1), INTENT(INOUT):: yq2
41 
42! Output arguments
43!****************************************************************************************
44    REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefh
45    REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefm
46
47! Other local variables
48!****************************************************************************************
49    INTEGER                                    :: k, i, j
50    REAL, DIMENSION(klon,klev)                 :: ycoefm0, ycoefh0, yzlay, yteta
51    REAL, DIMENSION(klon,klev+1)               :: yzlev, q2diag, ykmm, ykmn, ykmq
52    REAL, DIMENSION(klon)                      :: yustar
53
54! Include
55!****************************************************************************************
56    INCLUDE "clesphys.h"
57    INCLUDE "indicesol.h"
58    INCLUDE "iniprint.h"
59    INCLUDE "compbl.h"
60    INCLUDE "YOETHF.h"
61    INCLUDE "YOMCST.h"
62
63
64!****************************************************************************************   
65! Calcul de coefficients de diffusion turbulent de l'atmosphere :
66! ycoefm(:,2:klev), ycoefh(:,2:klev)
67!
68!****************************************************************************************   
69
70    CALL coefkz(nsrf, knon, ypaprs, ypplay, &
71         ksta, ksta_ter, &
72         yts, yrugos, yu, yv, yt, yq, &
73         yqsurf, &
74         ycoefm, ycoefh)
75 
76!****************************************************************************************
77! Eventuelle recalcule des coeffeicients de diffusion turbulent de l'atmosphere :
78! ycoefm(:,2:klev), ycoefh(:,2:klev)
79!
80!****************************************************************************************
81
82    IF (iflag_pbl.EQ.1) THEN
83       CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, &
84            ycoefm0, ycoefh0)
85
86       DO k = 2, klev
87          DO i = 1, knon
88             ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
89             ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
90          ENDDO
91       ENDDO
92    ENDIF
93
94 
95!**************************************************************************************** 
96! Calcul d'une diffusion minimale pour les conditions tres stables
97!
98!****************************************************************************************
99    IF (ok_kzmin) THEN
100       CALL coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycdragm, &
101            ycoefm0,ycoefh0)
102       
103       DO k = 2, klev
104          DO i = 1, knon
105             ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
106             ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
107          ENDDO
108       ENDDO
109       
110    ENDIF
111
112 
113!****************************************************************************************
114! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
115!
116!****************************************************************************************
117
118    IF (iflag_pbl.GE.3) THEN
119
120       yzlay(1:knon,1)= &
121            RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1))) &
122            *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
123       DO k=2,klev
124          DO i = 1, knon
125             yzlay(i,k)= &
126                  yzlay(i,k-1)+RD*0.5*(yt(i,k-1)+yt(i,k)) &
127                  /ypaprs(i,k)*(ypplay(i,k-1)-ypplay(i,k))/RG
128          END DO
129       END DO
130
131       DO k=1,klev
132          DO i = 1, knon
133             yteta(i,k)= &
134                  yt(i,k)*(ypaprs(i,1)/ypplay(i,k))**RKAPPA &
135                  *(1.+0.61*yq(i,k))
136          END DO
137       END DO
138
139       yzlev(1:knon,1)=0.
140       yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
141       DO k=2,klev
142          DO i = 1, knon
143             yzlev(i,k)=0.5*(yzlay(i,k)+yzlay(i,k-1))
144          END DO
145       END DO
146
147!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
148!!$! Pour memoire, le papier Hourdin et al. 2002 a ete obtenur avec un
149!!$! bug sur les coefficients de surface :
150!!$!          ycdragh(1:knon) = ycoefm(1:knon,1)
151!!$!          ycdragm(1:knon) = ycoefh(1:knon,1)
152!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
153       CALL ustarhb(knon,yu,yv,ycdragm, yustar)
154     
155       IF (prt_level > 9) THEN
156          WRITE(lunout,*) 'USTAR = ',yustar
157       ENDIF
158         
159!   iflag_pbl peut etre utilise comme longuer de melange
160       IF (iflag_pbl.GE.18) THEN
161          CALL vdif_kcay(knon,dtime,RG,RD,ypaprs,yt, &
162               yzlev,yzlay,yu,yv,yteta, &
163               ycdragm,yq2,q2diag,ykmm,ykmn,yustar, &
164               iflag_pbl)
165       ELSE
166          CALL yamada4(knon,dtime,RG,RD,ypaprs,yt, &
167               yzlev,yzlay,yu,yv,yteta, &
168               ycdragm,yq2,ykmm,ykmn,ykmq,yustar, &
169               iflag_pbl)
170       ENDIF
171       
172       ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
173       ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
174               
175    ENDIF !(iflag_pbl.ge.3)
176
177  END SUBROUTINE coef_diff_turb
178!
179!****************************************************************************************
180!
181  SUBROUTINE coefkz(nsrf, knon, paprs, pplay, &
182       ksta, ksta_ter, &
183       ts, rugos, &
184       u,v,t,q, &
185       qsurf, &
186       pcfm, pcfh)
187   
188    USE dimphy
189 
190!======================================================================
191! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
192!           (une version strictement identique a l'ancien modele)
193! Objet: calculer le coefficient du frottement du sol (Cdrag) et les
194!        coefficients d'echange turbulent dans l'atmosphere.
195! Arguments:
196! nsrf-----input-I- indicateur de la nature du sol
197! knon-----input-I- nombre de points a traiter
198! paprs----input-R- pregssion a chaque intercouche (en Pa)
199! pplay----input-R- pression au milieu de chaque couche (en Pa)
200! ts-------input-R- temperature du sol (en Kelvin)
201! rugos----input-R- longeur de rugosite (en m)
202! u--------input-R- vitesse u
203! v--------input-R- vitesse v
204! t--------input-R- temperature (K)
205! q--------input-R- vapeur d'eau (kg/kg)
206!
207! pcfm-----output-R- coefficients a calculer (vitesse)
208! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
209!======================================================================
210    INCLUDE "YOETHF.h"
211    INCLUDE "FCTTRE.h"
212    INCLUDE "iniprint.h"
213    INCLUDE "indicesol.h"
214    INCLUDE "compbl.h"
215    INCLUDE "YOMCST.h"
216!
217! Arguments:
218!
219    INTEGER, INTENT(IN)                      :: knon, nsrf
220    REAL, INTENT(IN)                         :: ksta, ksta_ter
221    REAL, DIMENSION(klon), INTENT(IN)        :: ts
222    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
223    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
224    REAL, DIMENSION(klon,klev), INTENT(IN)   :: u, v, t, q
225    REAL, DIMENSION(klon), INTENT(IN)        :: rugos
226    REAL, DIMENSION(klon), INTENT(IN)        :: qsurf
227
228    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: pcfm, pcfh
229
230!
231! Local variables:
232!
233    INTEGER, DIMENSION(klon)    :: itop ! numero de couche du sommet de la couche limite
234!
235! Quelques constantes et options:
236!
237    REAL, PARAMETER :: cepdu2=0.1**2
238    REAL, PARAMETER :: CKAP=0.4
239    REAL, PARAMETER :: cb=5.0
240    REAL, PARAMETER :: cc=5.0
241    REAL, PARAMETER :: cd=5.0
242    REAL, PARAMETER :: clam=160.0
243    REAL, PARAMETER :: ratqs=0.05 ! largeur de distribution de vapeur d'eau
244    LOGICAL, PARAMETER :: richum=.TRUE. ! utilise le nombre de Richardson humide
245    REAL, PARAMETER :: ric=0.4 ! nombre de Richardson critique
246    REAL, PARAMETER :: prandtl=0.4
247    REAL kstable ! diffusion minimale (situation stable)
248    ! GKtest
249    ! PARAMETER (kstable=1.0e-10)
250!IM: 261103     REAL kstable_ter, kstable_sinon
251!IM: 211003 cf GK   PARAMETER (kstable_ter = 1.0e-6)
252!IM: 261103     PARAMETER (kstable_ter = 1.0e-8)
253!IM: 261103   PARAMETER (kstable_ter = 1.0e-10)
254!IM: 261103   PARAMETER (kstable_sinon = 1.0e-10)
255    ! fin GKtest
256    REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
257    INTEGER isommet ! le sommet de la couche limite
258    LOGICAL, PARAMETER :: tvirtu=.TRUE. ! calculer Ri d'une maniere plus performante
259    LOGICAL, PARAMETER :: opt_ec=.FALSE.! formule du Centre Europeen dans l'atmosphere
260
261!
262! Variables locales:
263    INTEGER i, k !IM 120704
264    REAL zgeop(klon,klev)
265    REAL zmgeom(klon)
266    REAL zri(klon)
267    REAL zl2(klon)
268    REAL zdphi, zdu2, ztvd, ztvu, zcdn
269    REAL zscf
270    REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
271    REAL z2geomf, zalh2, zalm2, zscfh, zscfm
272    REAL, PARAMETER :: t_coup=273.15
273    LOGICAL, PARAMETER :: check=.FALSE.
274!
275! contre-gradient pour la chaleur sensible: Kelvin/metre
276    REAL gamt(2:klev)
277
278    LOGICAL, SAVE :: appel1er=.TRUE.
279    !$OMP THREADPRIVATE(appel1er)
280!
281! Fonctions thermodynamiques et fonctions d'instabilite
282    REAL fsta, fins, x
283
284    fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
285    fins(x) = SQRT(1.0-18.0*x)
286
287    isommet=klev
288     
289    IF (appel1er) THEN
290       IF (prt_level > 9) THEN
291          WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
292          WRITE(lunout,*)'coefkz, richum:', richum
293          IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
294          WRITE(lunout,*)'coefkz, isommet:', isommet
295          WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
296          appel1er = .FALSE.
297       ENDIF
298    ENDIF
299!
300! Initialiser les sorties
301!
302    DO k = 1, klev
303       DO i = 1, knon
304          pcfm(i,k) = 0.0
305          pcfh(i,k) = 0.0
306       ENDDO
307    ENDDO
308    DO i = 1, knon
309       itop(i) = 0
310    ENDDO
311   
312!
313! Prescrire la valeur de contre-gradient
314!
315    IF (iflag_pbl.EQ.1) THEN
316       DO k = 3, klev
317          gamt(k) = -1.0E-03
318       ENDDO
319       gamt(2) = -2.5E-03
320    ELSE
321       DO k = 2, klev
322          gamt(k) = 0.0
323       ENDDO
324    ENDIF
325!IM cf JLD/ GKtest
326    IF ( nsrf .NE. is_oce ) THEN
327!IM 261103     kstable = kstable_ter
328       kstable = ksta_ter
329    ELSE
330!IM 261103     kstable = kstable_sinon
331       kstable = ksta
332    ENDIF
333!IM cf JLD/ GKtest fin
334
335!
336! Calculer les geopotentiels de chaque couche
337!
338    DO i = 1, knon
339       zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) &
340            * (paprs(i,1)-pplay(i,1))
341    ENDDO
342    DO k = 2, klev
343       DO i = 1, knon
344          zgeop(i,k) = zgeop(i,k-1) &
345               + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) &
346               * (pplay(i,k-1)-pplay(i,k))
347       ENDDO
348    ENDDO
349
350!
351! Calculer les coefficients turbulents dans l'atmosphere
352!
353    DO i = 1, knon
354       itop(i) = isommet
355    ENDDO
356
357
358    DO k = 2, isommet
359       DO i = 1, knon
360          zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 &
361               +(v(i,k)-v(i,k-1))**2)
362          zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
363          zdphi =zmgeom(i) / 2.0
364          zt = (t(i,k)+t(i,k-1)) * 0.5
365          zq = (q(i,k)+q(i,k-1)) * 0.5
366
367!
368! Calculer Qs et dQs/dT:
369!
370          IF (thermcep) THEN
371             zdelta = MAX(0.,SIGN(1.,RTT-zt))
372             zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta) &
373                  + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
374             zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
375             zqs = MIN(0.5,zqs)
376             zcor = 1./(1.-RETV*zqs)
377             zqs = zqs*zcor
378             zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
379          ELSE
380             IF (zt .LT. t_coup) THEN
381                zqs = qsats(zt) / pplay(i,k)
382                zdqs = dqsats(zt,zqs)
383             ELSE
384                zqs = qsatl(zt) / pplay(i,k)
385                zdqs = dqsatl(zt,zqs)
386             ENDIF
387          ENDIF
388!
389!           calculer la fraction nuageuse (processus humide):
390!
391          if (zq /= 0.) then
392             zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
393          else
394             zfr = 0.
395          end if
396          zfr = MAX(0.0,MIN(1.0,zfr))
397          IF (.NOT.richum) zfr = 0.0
398!
399!           calculer le nombre de Richardson:
400!
401          IF (tvirtu) THEN
402             ztvd =( t(i,k) &
403                  + zdphi/RCPD/(1.+RVTMP2*zq) &
404                  *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
405                  )*(1.+RETV*q(i,k))
406             ztvu =( t(i,k-1) &
407                  - zdphi/RCPD/(1.+RVTMP2*zq) &
408                  *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
409                  )*(1.+RETV*q(i,k-1))
410             zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
411             zri(i) = zri(i) &
412                  + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
413                  *(paprs(i,k)/101325.0)**RKAPPA &
414                  /(zdu2*0.5*(ztvd+ztvu))
415
416          ELSE ! calcul de Ridchardson compatible LMD5
417
418             zri(i) =(RCPD*(t(i,k)-t(i,k-1)) &
419                  -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k) &
420                  *(pplay(i,k)-pplay(i,k-1)) &
421                  )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))
422             zri(i) = zri(i) + &
423                  zmgeom(i)*zmgeom(i)*gamt(k)/RG &
424                  *(paprs(i,k)/101325.0)**RKAPPA &
425                  /(zdu2*0.5*(t(i,k-1)+t(i,k)))
426          ENDIF
427!
428!           finalement, les coefficients d'echange sont obtenus:
429!
430          zcdn=SQRT(zdu2) / zmgeom(i) * RG
431
432          IF (opt_ec) THEN
433             z2geomf=zgeop(i,k-1)+zgeop(i,k)
434             zalm2=(0.5*ckap/RG*z2geomf &
435                  /(1.+0.5*ckap/rg/clam*z2geomf))**2
436             zalh2=(0.5*ckap/rg*z2geomf &
437                  /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
438             IF (zri(i).LT.0.0) THEN  ! situation instable
439                zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3 &
440                     / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
441                zscf = SQRT(-zri(i)*zscf)
442                zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
443                zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
444                pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
445                pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
446             ELSE ! situation stable
447                zscf=SQRT(1.+cd*zri(i))
448                pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
449                pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
450             ENDIF
451          ELSE
452             zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) &
453                  /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
454             pcfm(i,k)=SQRT(MAX(zcdn*zcdn*(ric-zri(i))/ric, kstable))
455             pcfm(i,k)= zl2(i)* pcfm(i,k)
456             pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
457          ENDIF
458       ENDDO
459    ENDDO
460
461!
462! Au-dela du sommet, pas de diffusion turbulente:
463!
464    DO i = 1, knon
465       IF (itop(i)+1 .LE. klev) THEN
466          DO k = itop(i)+1, klev
467             pcfh(i,k) = 0.0
468             pcfm(i,k) = 0.0
469          ENDDO
470       ENDIF
471    ENDDO
472     
473  END SUBROUTINE coefkz
474!
475!****************************************************************************************
476!
477  SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, &
478       pcfm, pcfh)
479
480    USE dimphy
481
482!======================================================================
483! J'introduit un peu de diffusion sauf dans les endroits
484! ou une forte inversion est presente
485! On peut dire qu'il represente la convection peu profonde
486!
487! Arguments:
488! nsrf-----input-I- indicateur de la nature du sol
489! knon-----input-I- nombre de points a traiter
490! paprs----input-R- pression a chaque intercouche (en Pa)
491! pplay----input-R- pression au milieu de chaque couche (en Pa)
492! t--------input-R- temperature (K)
493!
494! pcfm-----output-R- coefficients a calculer (vitesse)
495! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
496!======================================================================
497!
498! Arguments:
499!
500    INTEGER, INTENT(IN)                       :: knon, nsrf
501    REAL, DIMENSION(klon, klev+1), INTENT(IN) ::  paprs
502    REAL, DIMENSION(klon, klev), INTENT(IN)   ::  pplay
503    REAL, DIMENSION(klon, klev), INTENT(IN)   :: t(klon,klev)
504   
505    REAL, DIMENSION(klon, klev), INTENT(OUT)  :: pcfm, pcfh
506!
507! Quelques constantes et options:
508!
509    REAL, PARAMETER :: prandtl=0.4
510    REAL, PARAMETER :: kstable=0.002
511!   REAL, PARAMETER :: kstable=0.001
512    REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
513    REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
514!    PARAMETER (seuil=-0.04)
515!    PARAMETER (seuil=-0.06)
516!    PARAMETER (seuil=-0.09)
517
518!
519! Variables locales:
520!
521    INTEGER i, k, invb(knon)
522    REAL zl2(knon)
523    REAL zdthmin(knon), zdthdp
524
525    INCLUDE "indicesol.h"
526    INCLUDE "YOMCST.h"
527!
528! Initialiser les sorties
529!
530    DO k = 1, klev
531       DO i = 1, knon
532          pcfm(i,k) = 0.0
533          pcfh(i,k) = 0.0
534       ENDDO
535    ENDDO
536
537!
538! Chercher la zone d'inversion forte
539!
540    DO i = 1, knon
541       invb(i) = klev
542       zdthmin(i)=0.0
543    ENDDO
544    DO k = 2, klev/2-1
545       DO i = 1, knon
546          zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) &
547               - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
548          zdthdp = zdthdp * 100.0
549          IF (pplay(i,k).GT.0.8*paprs(i,1) .AND. &
550               zdthdp.LT.zdthmin(i) ) THEN
551             zdthmin(i) = zdthdp
552             invb(i) = k
553          ENDIF
554       ENDDO
555    ENDDO
556
557!
558! Introduire une diffusion:
559!
560    IF ( nsrf.EQ.is_oce ) THEN
561       DO k = 2, klev
562          DO i = 1, knon
563!IM cf FH/GK   IF ( (nsrf.NE.is_oce) .OR.  ! si ce n'est pas sur l'ocean
564!IM cf FH/GK  .     (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
565      !IM cf JLD/ GKtest TERkz2
566      ! IF ( (nsrf.EQ.is_ter) .OR.  ! si on est sur la terre
567      ! fin GKtest
568
569
570! s'il n'y a pas d'inversion ou si l'inversion est trop faible
571!          IF ( (nsrf.EQ.is_oce) .AND. &
572             IF ( (invb(i).EQ.klev) .OR. (zdthmin(i).GT.seuil) ) THEN
573                zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1)) &
574                     /(paprs(i,2)-paprs(i,klev+1)) ))**2
575                pcfm(i,k)= zl2(i)* kstable
576                pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
577             ENDIF
578          ENDDO
579       ENDDO
580    ENDIF
581
582  END SUBROUTINE coefkz2
583!
584!****************************************************************************************
585!
586END MODULE coef_diff_turb_mod
Note: See TracBrowser for help on using the repository browser.