My Project
 All Classes Files Functions Variables Macros
coef_diff_turb_mod.F90
Go to the documentation of this file.
1 !
2 MODULE 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 
10 CONTAINS
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.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 
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  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 
482 !======================================================================
483 ! J'introduit un peu de diffusion sauf dans les endroits
484 ! ou une forte inversion est presente
485 ! On peut dire qu'il represente la convection peu profonde
486 !
487 ! Arguments:
488 ! nsrf-----input-I- indicateur de la nature du sol
489 ! knon-----input-I- nombre de points a traiter
490 ! paprs----input-R- pression a chaque intercouche (en Pa)
491 ! pplay----input-R- pression au milieu de chaque couche (en Pa)
492 ! t--------input-R- temperature (K)
493 !
494 ! pcfm-----output-R- coefficients a calculer (vitesse)
495 ! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
496 !======================================================================
497 !
498 ! Arguments:
499 !
500  INTEGER, INTENT(IN) :: knon, nsrf
501  REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs
502  REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
503  REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon,klev)
504 
505  REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh
506 !
507 ! Quelques constantes et options:
508 !
509  REAL, PARAMETER :: prandtl=0.4
510  REAL, PARAMETER :: kstable=0.002
511 ! REAL, PARAMETER :: kstable=0.001
512  REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
513  REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
514 ! PARAMETER (seuil=-0.04)
515 ! PARAMETER (seuil=-0.06)
516 ! PARAMETER (seuil=-0.09)
517 
518 !
519 ! Variables locales:
520 !
521  INTEGER i, k, invb(knon)
522  REAL zl2(knon)
523  REAL zdthmin(knon), zdthdp
524 
525  include "indicesol.h"
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 !
586 END MODULE coef_diff_turb_mod