source: lmdz_wrf/WRFV3/lmdz/coef_diff_turb_mod.F90 @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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