source: LMDZ6/trunk/libf/phylmd/yamada_c.F90 @ 5451

Last change on this file since 5451 was 5451, checked in by fhourdin, 35 hours ago

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