source: LMDZ5/branches/LF-private/libf/phylmd/coef_diff_turb_mod.F90 @ 5425

Last change on this file since 5425 was 1785, checked in by Ehouarn Millour, 11 years ago

Transformation de l'include indicesol.h en un module indice_sol_mod et modification des appels dans tous les fichiers concernés.
Aucun changement des résultats ni des sorties du modèle vs 1784.
UG

...................................................

Replacement of the indicesol.h include by a module named indice_sol_mod. Modification of the calls in every affected files.
Results and outputs of simulations are unchanged in comparison with rev 1784.
UG

  • 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    USE indice_sol_mod
20!
21! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
22! atmosphere
23! NB! No values are calculated between surface and the first model layer.
24!     ycoefm(:,1) and ycoefh(:,1) are not valid !!!
25!
26!
27! Input arguments
28!****************************************************************************************
29    REAL, INTENT(IN)                           :: dtime
30    INTEGER, INTENT(IN)                        :: nsrf, knon
31    INTEGER, DIMENSION(klon), INTENT(IN)       :: ni
32    REAL, DIMENSION(klon,klev+1), INTENT(IN)   :: ypaprs
33    REAL, DIMENSION(klon,klev), INTENT(IN)     :: ypplay
34    REAL, DIMENSION(klon,klev), INTENT(IN)     :: yu, yv
35    REAL, DIMENSION(klon,klev), INTENT(IN)     :: yq, yt
36    REAL, DIMENSION(klon), INTENT(IN)          :: yts, yrugos, yqsurf
37    REAL, DIMENSION(klon), INTENT(IN)          :: ycdragm
38
39! InOutput arguments
40!****************************************************************************************
41    REAL, DIMENSION(klon,klev+1), INTENT(INOUT):: yq2
42 
43! Output arguments
44!****************************************************************************************
45    REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefh
46    REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefm
47
48! Other local variables
49!****************************************************************************************
50    INTEGER                                    :: k, i, j
51    REAL, DIMENSION(klon,klev)                 :: ycoefm0, ycoefh0, yzlay, yteta
52    REAL, DIMENSION(klon,klev+1)               :: yzlev, q2diag, ykmm, ykmn, ykmq
53    REAL, DIMENSION(klon)                      :: yustar
54
55! Include
56!****************************************************************************************
57    INCLUDE "clesphys.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.31) 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 IF (iflag_pbl<20) THEN
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    USE indice_sol_mod
190 
191!======================================================================
192! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
193!           (une version strictement identique a l'ancien modele)
194! Objet: calculer le coefficient du frottement du sol (Cdrag) et les
195!        coefficients d'echange turbulent dans l'atmosphere.
196! Arguments:
197! nsrf-----input-I- indicateur de la nature du sol
198! knon-----input-I- nombre de points a traiter
199! paprs----input-R- pregssion a chaque intercouche (en Pa)
200! pplay----input-R- pression au milieu de chaque couche (en Pa)
201! ts-------input-R- temperature du sol (en Kelvin)
202! rugos----input-R- longeur de rugosite (en m)
203! u--------input-R- vitesse u
204! v--------input-R- vitesse v
205! t--------input-R- temperature (K)
206! q--------input-R- vapeur d'eau (kg/kg)
207!
208! pcfm-----output-R- coefficients a calculer (vitesse)
209! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
210!======================================================================
211    INCLUDE "YOETHF.h"
212    INCLUDE "FCTTRE.h"
213    INCLUDE "iniprint.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    USE indice_sol_mod
482
483!======================================================================
484! J'introduit un peu de diffusion sauf dans les endroits
485! ou une forte inversion est presente
486! On peut dire qu'il represente la convection peu profonde
487!
488! Arguments:
489! nsrf-----input-I- indicateur de la nature du sol
490! knon-----input-I- nombre de points a traiter
491! paprs----input-R- pression a chaque intercouche (en Pa)
492! pplay----input-R- pression au milieu de chaque couche (en Pa)
493! t--------input-R- temperature (K)
494!
495! pcfm-----output-R- coefficients a calculer (vitesse)
496! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
497!======================================================================
498!
499! Arguments:
500!
501    INTEGER, INTENT(IN)                       :: knon, nsrf
502    REAL, DIMENSION(klon, klev+1), INTENT(IN) ::  paprs
503    REAL, DIMENSION(klon, klev), INTENT(IN)   ::  pplay
504    REAL, DIMENSION(klon, klev), INTENT(IN)   :: t(klon,klev)
505   
506    REAL, DIMENSION(klon, klev), INTENT(OUT)  :: pcfm, pcfh
507!
508! Quelques constantes et options:
509!
510    REAL, PARAMETER :: prandtl=0.4
511    REAL, PARAMETER :: kstable=0.002
512!   REAL, PARAMETER :: kstable=0.001
513    REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
514    REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
515!    PARAMETER (seuil=-0.04)
516!    PARAMETER (seuil=-0.06)
517!    PARAMETER (seuil=-0.09)
518
519!
520! Variables locales:
521!
522    INTEGER i, k, invb(knon)
523    REAL zl2(knon)
524    REAL zdthmin(knon), zdthdp
525
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.