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

Last change on this file since 4192 was 2952, checked in by Laurent Fairhead, 7 years ago

Parametrization of drag by copses
Need version 4465 of ORCHIDEE at least

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