source: LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90 @ 3982

Last change on this file since 3982 was 3577, checked in by fhourdin, 5 years ago

Corrections d'erreurs identifiees en debug ancienne physique sur Jean-Zay
Vive le poihl interactif.

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