source: LMDZ6/branches/Amaury_dev/libf/phylmd/yamada_c.F90 @ 5116

Last change on this file since 5116 was 5116, checked in by abarral, 2 months ago

rename modules properly lmdz_*
move ismin, ismax, minmax into new lmdz_libmath.f90
(lint) uppercase fortran keywords
  • 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
File size: 15.4 KB
Line 
1
2! $Header$
3
4      SUBROUTINE yamada_c(ngrid,timestep,plev,play &
5     ,pu,pv,pt,d_u,d_v,d_t,cd,q2,km,kn,kq,d_t_diss,ustar &
6     ,iflag_pbl)
7      USE dimphy, ONLY: klon, klev
8      USE lmdz_print_control, ONLY: prt_level
9      USE lmdz_ioipsl_getin_p, ONLY: getin_p
10
11      IMPLICIT NONE
12      INCLUDE "YOMCST.h"
13
14! timestep : pas de temps
15! g  : g
16! zlev : altitude a chaque niveau (interface inferieure de la couche
17!        de meme indice)
18! zlay : altitude au centre de chaque couche
19! u,v : vitesse au centre de chaque couche
20!       (en entree : la valeur au debut du pas de temps)
21! teta : temperature potentielle au centre de chaque couche
22!        (en entree : la valeur au debut du pas de temps)
23! cd : cdrag
24!      (en entree : la valeur au debut du pas de temps)
25! q2 : $q^2$ au bas de chaque couche
26!      (en entree : la valeur au debut du pas de temps)
27!      (en sortie : la valeur a la fin du pas de temps)
28! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
29!      couche)
30!      (en sortie : la valeur a la fin du pas de temps)
31! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
32!      (en sortie : la valeur a la fin du pas de temps)
33
34!  iflag_pbl doit valoir entre 6 et 9
35!      l=6, on prend  systematiquement une longueur d'equilibre
36!    iflag_pbl=6 : MY 2.0
37!    iflag_pbl=7 : MY 2.0.Fournier
38!    iflag_pbl=8/9 : MY 2.5
39!       iflag_pbl=8 with special obsolete treatments for convergence
40!       with Cmpi5 NPv3.1 simulations
41!    iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact
42!    iflag_pbl=12 = 11 with vertical diffusion off q2
43
44!  2013/04/01 (FH hourdin@lmd.jussieu.fr)
45!     Correction for very stable PBLs (iflag_pbl=10 and 11)
46!     iflag_pbl=8 converges numerically with NPv3.1
47!     iflag_pbl=11 -> the model starts with NP from start files created by ce0l
48!                  -> the model can run with longer time-steps.
49!.......................................................................
50
51      REAL, DIMENSION(klon,klev) :: d_u,d_v,d_t
52      REAL, DIMENSION(klon,klev) :: pu,pv,pt
53      REAL, DIMENSION(klon,klev) :: d_t_diss
54
55      REAL timestep
56      real plev(klon,klev+1)
57      real play(klon,klev)
58      real ustar(klon)
59      real kmin,qmin,pblhmin(klon),coriol(klon)
60      REAL zlev(klon,klev+1)
61      REAL zlay(klon,klev)
62      REAL zu(klon,klev)
63      REAL zv(klon,klev)
64      REAL zt(klon,klev)
65      REAL teta(klon,klev)
66      REAL cd(klon)
67      REAL q2(klon,klev+1),qpre
68      REAL unsdz(klon,klev)
69      REAL unsdzdec(klon,klev+1)
70
71      REAL km(klon,klev)
72      REAL kmpre(klon,klev+1),tmp2
73      REAL mpre(klon,klev+1)
74      REAL kn(klon,klev)
75      REAL kq(klon,klev)
76      real ff(klon,klev+1),delta(klon,klev+1)
77      real aa(klon,klev+1),aa0,aa1
78      integer iflag_pbl,ngrid
79      integer nlay,nlev
80
81      logical first
82      integer ipas
83      save first,ipas
84!FH/IM     data first,ipas/.TRUE.,0/
85      data first,ipas/.FALSE.,0/
86!$OMP THREADPRIVATE( first,ipas)
87       INTEGER, SAVE :: iflag_tke_diff=0
88!$OMP THREADPRIVATE(iflag_tke_diff)
89
90
91      integer ig,k
92
93
94      real ri,zrif,zalpha,zsm,zsn
95      real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
96
97      real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
98      REAL, DIMENSION(klon,klev+1) :: km2,kn2,sqrtq
99      real dtetadz(klon,klev+1)
100      real m2cstat,mcstat,kmcstat
101      real l(klon,klev+1)
102      real leff(klon,klev+1)
103      real,allocatable,save :: l0(:)
104!$OMP THREADPRIVATE(l0)     
105      real sq(klon),sqz(klon),zz(klon,klev+1)
106      integer iter
107
108      real ric,rifc,b1,kap
109      save ric,rifc,b1,kap
110      data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
111!$OMP THREADPRIVATE(ric,rifc,b1,kap)
112      real frif,falpha,fsm
113      real fl,zzz,zl0,zq2,zn2
114
115      real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
116      real lyam(klon,klev),knyam(klon,klev)
117      real w2yam(klon,klev),t2yam(klon,klev)
118      logical,save :: firstcall=.TRUE.
119!$OMP THREADPRIVATE(firstcall)       
120      CHARACTER(len=20),PARAMETER :: modname="yamada_c"
121REAL, DIMENSION(klon,klev+1) :: fluxu,fluxv,fluxt
122REAL, DIMENSION(klon,klev+1) :: dddu,dddv,dddt
123REAL, DIMENSION(klon,klev) :: exner,masse
124REAL, DIMENSION(klon,klev+1) :: masseb,q2old,q2neg
125      LOGICAL okiophys
126
127      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
128      falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
129      fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
130      fl(zzz,zl0,zq2,zn2)= &
131       max(min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig)) &
132       ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
133
134
135      okiophys=klon==1
136      if (firstcall) THEN
137        CALL getin_p('iflag_tke_diff',iflag_tke_diff)
138        allocate(l0(klon))
139#define IOPHYS
140#ifdef IOPHYS
141!        CALL iophys_ini(timestep)
142#endif
143        firstcall=.FALSE.
144      endif
145
146   IF (ngrid<=0) RETURN ! Bizarre : on n a pas ce probeleme pour coef_diff_turb
147
148#ifdef IOPHYS
149IF (okiophys) THEN
150CALL iophys_ecrit('q2i',klev,'q2 debut my','m2/s2',q2(:,1:klev))
151CALL iophys_ecrit('kmi',klev,'Kz debut my','m/s2',km(:,1:klev))
152END IF
153#endif
154
155      nlay=klev
156      nlev=klev+1
157
158
159!-------------------------------------------------------------------------
160! Computation of conservative source terms from the turbulent tendencies
161!-------------------------------------------------------------------------
162
163
164   zalpha=0.5 ! Anciennement 0.5. Essayer de voir pourquoi ?
165   zu(:,:)=pu(:,:)+zalpha*d_u(:,:)
166   zv(:,:)=pv(:,:)+zalpha*d_v(:,:)
167   zt(:,:)=pt(:,:)+zalpha*d_t(:,:)
168
169   do k=1,klev
170      exner(:,k)=(play(:,k)/plev(:,1))**RKAPPA
171      masse(:,k)=(plev(:,k)-plev(:,k+1))/RG
172      teta(:,k)=zt(:,k)/exner(:,k)
173   enddo
174
175! Atmospheric mass at layer interfaces, where the TKE is computed
176   masseb(:,:)=0.
177   do k=1,klev
178      masseb(:,k)=masseb(:,k)+masse(:,k)
179      masseb(:,k+1)=masseb(:,k+1)+masse(:,k)
180    enddo
181    masseb(:,:)=0.5*masseb(:,:)
182
183   zlev(:,1)=0.
184   zlay(:,1)=RCPD*teta(:,1)*(1.-exner(:,1))
185   do k=1,klev-1
186      zlay(:,k+1)=zlay(:,k)+0.5*RCPD*(teta(:,k)+teta(:,k+1))*(exner(:,k)-exner(:,k+1))/RG
187      zlev(:,k)=0.5*(zlay(:,k)+zlay(:,k+1)) ! PASBO
188   enddo
189
190   fluxu(:,klev+1)=0.
191   fluxv(:,klev+1)=0.
192   fluxt(:,klev+1)=0.
193
194   do k=klev,1,-1
195      fluxu(:,k)=fluxu(:,k+1)+masse(:,k)*d_u(:,k)
196      fluxv(:,k)=fluxv(:,k+1)+masse(:,k)*d_v(:,k)
197      fluxt(:,k)=fluxt(:,k+1)+masse(:,k)*d_t(:,k)/exner(:,k) ! Flux de theta
198   enddo
199
200   dddu(:,1)=2*zu(:,1)*fluxu(:,1)
201   dddv(:,1)=2*zv(:,1)*fluxv(:,1)
202   dddt(:,1)=(exner(:,1)-1.)*fluxt(:,1)
203
204   do k=2,klev
205      dddu(:,k)=(zu(:,k)-zu(:,k-1))*fluxu(:,k)
206      dddv(:,k)=(zv(:,k)-zv(:,k-1))*fluxv(:,k)
207      dddt(:,k)=(exner(:,k)-exner(:,k-1))*fluxt(:,k)
208   enddo
209   dddu(:,klev+1)=0.
210   dddv(:,klev+1)=0.
211   dddt(:,klev+1)=0.
212
213#ifdef IOPHYS
214IF (okiophys) THEN
215      CALL iophys_ecrit('zlay',klev,'Geop','m',zlay)
216      CALL iophys_ecrit('teta',klev,'teta','K',teta)
217      CALL iophys_ecrit('temp',klev,'temp','K',zt)
218      CALL iophys_ecrit('pt',klev,'temp','K',pt)
219      CALL iophys_ecrit('pu',klev,'u','m/s',pu)
220      CALL iophys_ecrit('pv',klev,'v','m/s',pv)
221      CALL iophys_ecrit('d_u',klev,'d_u','m/s2',d_u)
222      CALL iophys_ecrit('d_v',klev,'d_v','m/s2',d_v)
223      CALL iophys_ecrit('d_t',klev,'d_t','K/s',d_t)
224      CALL iophys_ecrit('exner',klev,'exner','',exner)
225      CALL iophys_ecrit('masse',klev,'masse','',masse)
226      CALL iophys_ecrit('masseb',klev,'masseb','',masseb)
227END IF
228#endif
229
230
231
232      ipas=ipas+1
233
234
235!.......................................................................
236!  les increments verticaux
237!.......................................................................
238
239!!!!!! allerte !!!!!c
240!!!!!! zlev n'est pas declare a nlev !!!!!c
241!!!!!! ---->
242                                                      DO ig=1,ngrid
243            zlev(ig,nlev)=zlay(ig,nlay) &
244               +( zlay(ig,nlay) - zlev(ig,nlev-1) )
245                                                      ENDDO
246!!!!!! <----
247!!!!!! allerte !!!!!c
248
249      DO k=1,nlay
250                                                      DO ig=1,ngrid
251        unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
252                                                      ENDDO
253      ENDDO
254                                                      DO ig=1,ngrid
255      unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
256                                                      ENDDO
257      DO k=2,nlay
258                                                      DO ig=1,ngrid
259        unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
260                                                     ENDDO
261      ENDDO
262                                                      DO ig=1,ngrid
263      unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
264                                                     ENDDO
265
266!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
267! Computing M^2, N^2, Richardson numbers, stability functions
268!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
269
270
271      do k=2,klev
272                                                          do ig=1,ngrid
273         dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
274         m2(ig,k)=((zu(ig,k)-zu(ig,k-1))**2+(zv(ig,k)-zv(ig,k-1))**2)/(dz(ig,k)*dz(ig,k))
275         dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
276         n2(ig,k)=RG*2.*dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
277!        n2(ig,k)=0.
278         ri=n2(ig,k)/max(m2(ig,k),1.e-10)
279         if (ri<ric) THEN
280            rif(ig,k)=frif(ri)
281         else
282            rif(ig,k)=rifc
283         endif
284         IF(rif(ig,k)<0.16) THEN
285            alpha(ig,k)=falpha(rif(ig,k))
286            sm(ig,k)=fsm(rif(ig,k))
287         else
288            alpha(ig,k)=1.12
289            sm(ig,k)=0.085
290         endif
291         zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
292                                                          enddo
293      enddo
294
295
296
297!====================================================================
298!  Computing the mixing length
299!====================================================================
300
301!   Mise a jour de l0
302      if (iflag_pbl==8.or.iflag_pbl==10) THEN
303!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
304! Iterative computation of l0
305! This version is kept for iflag_pbl only for convergence
306! with NPv3.1 Cmip5 simulations
307!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
308
309                                                          do ig=1,ngrid
310      sq(ig)=1.e-10
311      sqz(ig)=1.e-10
312                                                          enddo
313      do k=2,klev-1
314                                                          do ig=1,ngrid
315        zq=sqrt(q2(ig,k))
316        sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
317        sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
318                                                          enddo
319      enddo
320                                                          do ig=1,ngrid
321      l0(ig)=0.2*sqz(ig)/sq(ig)
322                                                          enddo
323      do k=2,klev
324                                                          do ig=1,ngrid
325         l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
326                                                          enddo
327      enddo
328!     PRINT*,'L0 cas 8 ou 10 ',l0
329
330      else
331
332!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
333! In all other case, the assymptotic mixing length l0 is imposed (100m)
334!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
335
336          l0(:)=150.
337          do k=2,klev
338                                                          do ig=1,ngrid
339             l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
340                                                          enddo
341          enddo
342!     PRINT*,'L0 cas autres ',l0
343
344      endif
345
346
347#ifdef IOPHYS
348IF (okiophys) THEN
349CALL iophys_ecrit('rif',klev,'Flux Richardson','m',rif(:,1:klev))
350CALL iophys_ecrit('m2',klev,'m2 ','m/s',m2(:,1:klev))
351CALL iophys_ecrit('Km2app',klev,'m2 conserv','m/s',km(:,1:klev)*m2(:,1:klev))
352CALL iophys_ecrit('Km',klev,'Km','m2/s',km(:,1:klev))
353END IF
354#endif
355
356
357IF (iflag_pbl<20) THEN
358      ! For diagnostics only
359      RETURN
360
361ELSE
362
363!  PRINT*,'OK1'
364
365! Evolution of TKE under source terms K M2 and K N2
366   leff(:,:)=max(l(:,:),1.)
367
368!##################################################################
369!#  IF (iflag_pbl==29) THEN
370!#     stop 'Ne pas utiliser iflag_pbl=29'
371!#     km2(:,:)=km(:,:)*m2(:,:)
372!#     kn2(:,:)=kn2(:,:)*rif(:,:)
373!#  ELSEIF (iflag_pbl==25) THEN
374! VERSION AVEC LA TKE EN MILIEU DE COUCHE
375!#     stop 'Ne pas utiliser iflag_pbl=25'
376!#     DO k=1,klev
377!#        km2(:,k)=-0.5*(dddu(:,k)+dddv(:,k)+dddu(:,k+1)+dddv(:,k+1)) &
378!#        &        /(masse(:,k)*timestep)
379!#        kn2(:,k)=rcpd*0.5*(dddt(:,k)+dddt(:,k+1))/(masse(:,k)*timestep)
380!#        leff(:,k)=0.5*(leff(:,k)+leff(:,k+1))
381!#     ENDDO
382!#     km2(:,klev+1)=0. ; kn2(:,klev+1)=0.
383!#  ELSE
384!#################################################################
385
386      km2(:,:)=-(dddu(:,:)+dddv(:,:))/(masseb(:,:)*timestep)
387      kn2(:,:)=rcpd*dddt(:,:)/(masseb(:,:)*timestep)
388!   ENDIF
389   q2neg(:,:)=q2(:,:)+timestep*(km2(:,:)-kn2(:,:))
390   q2(:,:)=min(max(q2neg(:,:),1.e-10),1.e4)
391
392 
393#ifdef IOPHYS
394IF (okiophys) THEN
395      CALL iophys_ecrit('km2',klev,'m2 conserv','m/s',km2(:,1:klev))
396      CALL iophys_ecrit('kn2',klev,'n2 conserv','m/s',kn2(:,1:klev))
397END IF
398#endif
399
400! Dissipation of TKE
401   q2old(:,:)=q2(:,:)
402   q2(:,:)=1./(1./sqrt(q2(:,:))+timestep/(2*leff(:,:)*b1))
403   q2(:,:)=q2(:,:)*q2(:,:)
404!  IF (iflag_pbl<=24) THEN
405      DO k=1,klev
406         d_t_diss(:,k)=(masseb(:,k)*(q2neg(:,k)-q2(:,k))+masseb(:,k+1)*(q2neg(:,k+1)-q2(:,k+1)))/(2.*rcpd*masse(:,k))
407      ENDDO
408
409!###################################################################
410!  ELSE IF (iflag_pbl<=27) THEN
411!     DO k=1,klev
412!        d_t_diss(:,k)=(q2neg(:,k)-q2(:,k))/rcpd
413!     ENDDO
414!  ENDIF
415!  PRINT*,'iflag_pbl ',d_t_diss
416!###################################################################
417
418
419! Compuation of stability functions
420!   IF (iflag_pbl/=29) THEN
421      DO k=1,klev
422      DO ig=1,ngrid
423         IF (ABS(km2(ig,k))<=1.e-20) THEN
424            rif(ig,k)=0.
425         ELSE
426            rif(ig,k)=min(kn2(ig,k)/km2(ig,k),rifc)
427         ENDIF
428         IF (rif(ig,k)<0.16) THEN
429            alpha(ig,k)=falpha(rif(ig,k))
430            sm(ig,k)=fsm(rif(ig,k))
431         else
432            alpha(ig,k)=1.12
433            sm(ig,k)=0.085
434         endif
435      ENDDO
436      ENDDO
437!    ENDIF
438
439! Computation of turbulent diffusivities
440!  IF (25<=iflag_pbl.and.iflag_pbl<=28) THEN
441!    DO k=2,klev
442!       sqrtq(:,k)=sqrt(0.5*(q2(:,k)+q2(:,k-1)))
443!    ENDDO
444!  ELSE
445   kq(:,:)=0.
446   DO k=1,klev
447      ! Coefficient au milieu des couches pour diffuser la TKE
448      kq(:,k)=0.5*leff(:,k)*sqrt(q2(:,k))*0.2
449   ENDDO
450
451#ifdef IOPHYS
452IF (okiophys) THEN
453CALL iophys_ecrit('q2b',klev,'KTE inter','m2/s',q2(:,1:klev))
454END IF
455#endif
456
457  IF (iflag_tke_diff==1) THEN
458    CALL vdif_q2(timestep, RG, RD, ngrid, plev, pt, kq, q2)
459  ENDIF
460
461   km(:,:)=0.
462   kn(:,:)=0.
463   DO k=1,klev
464      km(:,k)=leff(:,k)*sqrt(q2(:,k))*sm(:,k)
465      kn(:,k)=km(:,k)*alpha(:,k)
466   ENDDO
467
468
469#ifdef IOPHYS
470IF (okiophys) THEN
471CALL iophys_ecrit('mixingl',klev,'Mixing length','m',leff(:,1:klev))
472CALL iophys_ecrit('rife',klev,'Flux Richardson','m',rif(:,1:klev))
473CALL iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
474CALL iophys_ecrit('q2neg',klev,'KTE non bornee','m2/s',q2neg(:,1:klev))
475CALL iophys_ecrit('alpha',klev,'alpha','',alpha(:,1:klev))
476CALL iophys_ecrit('sm',klev,'sm','',sm(:,1:klev))
477CALL iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
478CALL iophys_ecrit('kmf',klev,'Kz final','m2/s',km(:,1:klev))
479CALL iophys_ecrit('knf',klev,'Kz final','m2/s',kn(:,1:klev))
480CALL iophys_ecrit('kqf',klev,'Kz final','m2/s',kq(:,1:klev))
481END IF
482#endif
483
484
485ENDIF
486
487
488!  PRINT*,'OK2'
489      RETURN
490      END
Note: See TracBrowser for help on using the repository browser.