source: trunk/LMDZ.EARTH/libf/phylmd/coef_diff_turb_mod.F90 @ 135

Last change on this file since 135 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5

File size: 18.7 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.11) 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          zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
392          zfr = MAX(0.0,MIN(1.0,zfr))
393          IF (.NOT.richum) zfr = 0.0
394!
395!           calculer le nombre de Richardson:
396!
397          IF (tvirtu) THEN
398             ztvd =( t(i,k) &
399                  + zdphi/RCPD/(1.+RVTMP2*zq) &
400                  *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
401                  )*(1.+RETV*q(i,k))
402             ztvu =( t(i,k-1) &
403                  - zdphi/RCPD/(1.+RVTMP2*zq) &
404                  *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
405                  )*(1.+RETV*q(i,k-1))
406             zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
407             zri(i) = zri(i) &
408                  + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
409                  *(paprs(i,k)/101325.0)**RKAPPA &
410                  /(zdu2*0.5*(ztvd+ztvu))
411
412          ELSE ! calcul de Ridchardson compatible LMD5
413
414             zri(i) =(RCPD*(t(i,k)-t(i,k-1)) &
415                  -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k) &
416                  *(pplay(i,k)-pplay(i,k-1)) &
417                  )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))
418             zri(i) = zri(i) + &
419                  zmgeom(i)*zmgeom(i)*gamt(k)/RG &
420                  *(paprs(i,k)/101325.0)**RKAPPA &
421                  /(zdu2*0.5*(t(i,k-1)+t(i,k)))
422          ENDIF
423!
424!           finalement, les coefficients d'echange sont obtenus:
425!
426          zcdn=SQRT(zdu2) / zmgeom(i) * RG
427
428          IF (opt_ec) THEN
429             z2geomf=zgeop(i,k-1)+zgeop(i,k)
430             zalm2=(0.5*ckap/RG*z2geomf &
431                  /(1.+0.5*ckap/rg/clam*z2geomf))**2
432             zalh2=(0.5*ckap/rg*z2geomf &
433                  /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
434             IF (zri(i).LT.0.0) THEN  ! situation instable
435                zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3 &
436                     / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
437                zscf = SQRT(-zri(i)*zscf)
438                zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
439                zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
440                pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
441                pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
442             ELSE ! situation stable
443                zscf=SQRT(1.+cd*zri(i))
444                pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
445                pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
446             ENDIF
447          ELSE
448             zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) &
449                  /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
450             pcfm(i,k)=SQRT(MAX(zcdn*zcdn*(ric-zri(i))/ric, kstable))
451             pcfm(i,k)= zl2(i)* pcfm(i,k)
452             pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
453          ENDIF
454       ENDDO
455    ENDDO
456
457!
458! Au-dela du sommet, pas de diffusion turbulente:
459!
460    DO i = 1, knon
461       IF (itop(i)+1 .LE. klev) THEN
462          DO k = itop(i)+1, klev
463             pcfh(i,k) = 0.0
464             pcfm(i,k) = 0.0
465          ENDDO
466       ENDIF
467    ENDDO
468     
469  END SUBROUTINE coefkz
470!
471!****************************************************************************************
472!
473  SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, &
474       pcfm, pcfh)
475
476    USE dimphy
477
478!======================================================================
479! J'introduit un peu de diffusion sauf dans les endroits
480! ou une forte inversion est presente
481! On peut dire qu'il represente la convection peu profonde
482!
483! Arguments:
484! nsrf-----input-I- indicateur de la nature du sol
485! knon-----input-I- nombre de points a traiter
486! paprs----input-R- pression a chaque intercouche (en Pa)
487! pplay----input-R- pression au milieu de chaque couche (en Pa)
488! t--------input-R- temperature (K)
489!
490! pcfm-----output-R- coefficients a calculer (vitesse)
491! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
492!======================================================================
493!
494! Arguments:
495!
496    INTEGER, INTENT(IN)                       :: knon, nsrf
497    REAL, DIMENSION(klon, klev+1), INTENT(IN) ::  paprs
498    REAL, DIMENSION(klon, klev), INTENT(IN)   ::  pplay
499    REAL, DIMENSION(klon, klev), INTENT(IN)   :: t(klon,klev)
500   
501    REAL, DIMENSION(klon, klev), INTENT(OUT)  :: pcfm, pcfh
502!
503! Quelques constantes et options:
504!
505    REAL, PARAMETER :: prandtl=0.4
506    REAL, PARAMETER :: kstable=0.002
507!   REAL, PARAMETER :: kstable=0.001
508    REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
509    REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
510!    PARAMETER (seuil=-0.04)
511!    PARAMETER (seuil=-0.06)
512!    PARAMETER (seuil=-0.09)
513
514!
515! Variables locales:
516!
517    INTEGER i, k, invb(knon)
518    REAL zl2(knon)
519    REAL zdthmin(knon), zdthdp
520
521    INCLUDE "indicesol.h"
522    INCLUDE "YOMCST.h"
523!
524! Initialiser les sorties
525!
526    DO k = 1, klev
527       DO i = 1, knon
528          pcfm(i,k) = 0.0
529          pcfh(i,k) = 0.0
530       ENDDO
531    ENDDO
532
533!
534! Chercher la zone d'inversion forte
535!
536    DO i = 1, knon
537       invb(i) = klev
538       zdthmin(i)=0.0
539    ENDDO
540    DO k = 2, klev/2-1
541       DO i = 1, knon
542          zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) &
543               - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
544          zdthdp = zdthdp * 100.0
545          IF (pplay(i,k).GT.0.8*paprs(i,1) .AND. &
546               zdthdp.LT.zdthmin(i) ) THEN
547             zdthmin(i) = zdthdp
548             invb(i) = k
549          ENDIF
550       ENDDO
551    ENDDO
552
553!
554! Introduire une diffusion:
555!
556    IF ( nsrf.EQ.is_oce ) THEN
557       DO k = 2, klev
558          DO i = 1, knon
559!IM cf FH/GK   IF ( (nsrf.NE.is_oce) .OR.  ! si ce n'est pas sur l'ocean
560!IM cf FH/GK  .     (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
561      !IM cf JLD/ GKtest TERkz2
562      ! IF ( (nsrf.EQ.is_ter) .OR.  ! si on est sur la terre
563      ! fin GKtest
564
565
566! s'il n'y a pas d'inversion ou si l'inversion est trop faible
567!          IF ( (nsrf.EQ.is_oce) .AND. &
568             IF ( (invb(i).EQ.klev) .OR. (zdthmin(i).GT.seuil) ) THEN
569                zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1)) &
570                     /(paprs(i,2)-paprs(i,klev+1)) ))**2
571                pcfm(i,k)= zl2(i)* kstable
572                pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
573             ENDIF
574          ENDDO
575       ENDDO
576    ENDIF
577
578  END SUBROUTINE coefkz2
579!
580!****************************************************************************************
581!
582END MODULE coef_diff_turb_mod
Note: See TracBrowser for help on using the repository browser.