source: LMDZ4/trunk/libf/phylmd/coef_diff_turb_mod.F90 @ 862

Last change on this file since 862 was 793, checked in by Laurent Fairhead, 17 years ago

Modifications suite a la transformation des fichiers include pour
qu'ils soient compatibles a la fois au format fixe et au format libre
Un bon nombre de fichiers *.inc du coup disparaissent
LF

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