source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/yamada4.F @ 1299

Last change on this file since 1299 was 1299, checked in by Laurent Fairhead, 14 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.6 KB
Line 
1!
2! $Id: yamada4.F 1299 2010-01-20 14:27:21Z fairhead $
3!
4      SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp
5     s   ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
6     s   ,iflag_pbl)
7      use dimphy
8      IMPLICIT NONE
9#include "iniprint.h"
10c.......................................................................
11cym#include "dimensions.h"
12cym#include "dimphy.h"
13c.......................................................................
14c
15c dt : pas de temps
16c g  : g
17c zlev : altitude a chaque niveau (interface inferieure de la couche
18c        de meme indice)
19c zlay : altitude au centre de chaque couche
20c u,v : vitesse au centre de chaque couche
21c       (en entree : la valeur au debut du pas de temps)
22c teta : temperature potentielle au centre de chaque couche
23c        (en entree : la valeur au debut du pas de temps)
24c cd : cdrag
25c      (en entree : la valeur au debut du pas de temps)
26c q2 : $q^2$ au bas de chaque couche
27c      (en entree : la valeur au debut du pas de temps)
28c      (en sortie : la valeur a la fin du pas de temps)
29c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
30c      couche)
31c      (en sortie : la valeur a la fin du pas de temps)
32c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
33c      (en sortie : la valeur a la fin du pas de temps)
34c
35c  iflag_pbl doit valoir entre 6 et 9
36c      l=6, on prend  systematiquement une longueur d'equilibre
37c    iflag_pbl=6 : MY 2.0
38c    iflag_pbl=7 : MY 2.0.Fournier
39c    iflag_pbl=8 : MY 2.5
40c    iflag_pbl=9 : un test ?
41
42c.......................................................................
43      REAL dt,g,rconst
44      real plev(klon,klev+1),temp(klon,klev)
45      real ustar(klon)
46      real kmin,qmin,pblhmin(klon),coriol(klon)
47      REAL zlev(klon,klev+1)
48      REAL zlay(klon,klev)
49      REAL u(klon,klev)
50      REAL v(klon,klev)
51      REAL teta(klon,klev)
52      REAL cd(klon)
53      REAL q2(klon,klev+1),qpre
54      REAL unsdz(klon,klev)
55      REAL unsdzdec(klon,klev+1)
56
57      REAL km(klon,klev+1)
58      REAL kmpre(klon,klev+1),tmp2
59      REAL mpre(klon,klev+1)
60      REAL kn(klon,klev+1)
61      REAL kq(klon,klev+1)
62      real ff(klon,klev+1),delta(klon,klev+1)
63      real aa(klon,klev+1),aa0,aa1
64      integer iflag_pbl,ngrid
65
66
67      integer nlay,nlev
68cym      PARAMETER (nlay=klev)
69cym      PARAMETER (nlev=klev+1)
70
71      logical first
72      integer ipas
73      save first,ipas
74cFH/IM     data first,ipas/.true.,0/
75      data first,ipas/.false.,0/
76c$OMP THREADPRIVATE( first,ipas)
77
78      integer ig,k
79
80
81      real ri,zrif,zalpha,zsm,zsn
82      real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
83
84      real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
85      real dtetadz(klon,klev+1)
86      real m2cstat,mcstat,kmcstat
87      real l(klon,klev+1)
88      real,allocatable,save :: l0(:)
89c$OMP THREADPRIVATE(l0)     
90      real sq(klon),sqz(klon),zz(klon,klev+1)
91      integer iter
92
93      real ric,rifc,b1,kap
94      save ric,rifc,b1,kap
95      data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
96c$OMP THREADPRIVATE(ric,rifc,b1,kap)
97      real frif,falpha,fsm
98      real fl,zzz,zl0,zq2,zn2
99
100cym      real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
101cym     s  ,lyam(klon,klev),knyam(klon,klev)
102cym     s  ,w2yam(klon,klev),t2yam(klon,klev)
103      real,allocatable,save,dimension(:,:) :: rino,smyam,styam,lyam,
104     s                                        knyam,w2yam,t2yam
105cym      common/pbldiag/rino,smyam,styam,lyam,knyam,w2yam,t2yam
106c$OMP THREADPRIVATE(rino,smyam,styam,lyam,knyam,w2yam,t2yam)
107      logical,save :: firstcall=.true.
108
109      character (len=20) :: modname='yamada4'
110      character (len=80) :: abort_message
111
112c$OMP THREADPRIVATE(firstcall)       
113      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
114      falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
115      fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
116      fl(zzz,zl0,zq2,zn2)=
117     s     max(min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
118     s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
119
120
121      nlay=klev
122      nlev=klev+1
123     
124      if (firstcall) then
125        allocate(rino(klon,klev+1),smyam(klon,klev),styam(klon,klev))
126        allocate(lyam(klon,klev),knyam(klon,klev))
127        allocate(w2yam(klon,klev),t2yam(klon,klev))
128        allocate(l0(klon))
129        firstcall=.false.
130      endif
131
132
133      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
134              abort_message = 'probleme de coherence dans appel a MY'
135              CALL abort_gcm (modname,abort_message,1)
136      endif
137
138      ipas=ipas+1
139      if (0.eq.1.and.first) then
140      do ig=1,1000
141         ri=(ig-800.)/500.
142         if (ri.lt.ric) then
143            zrif=frif(ri)
144         else
145            zrif=rifc
146         endif
147         if(zrif.lt.0.16) then
148            zalpha=falpha(zrif)
149            zsm=fsm(zrif)
150         else
151            zalpha=1.12
152            zsm=0.085
153         endif
154c     print*,ri,rif,zalpha,zsm
155      enddo
156      endif
157
158c.......................................................................
159c  les increments verticaux
160c.......................................................................
161c
162c!!!!! allerte !!!!!c
163c!!!!! zlev n'est pas declare a nlev !!!!!c
164c!!!!! ---->
165                                                      DO ig=1,ngrid
166            zlev(ig,nlev)=zlay(ig,nlay)
167     &             +( zlay(ig,nlay) - zlev(ig,nlev-1) )
168                                                      ENDDO
169c!!!!! <----
170c!!!!! allerte !!!!!c
171c
172      DO k=1,nlay
173                                                      DO ig=1,ngrid
174        unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
175                                                      ENDDO
176      ENDDO
177                                                      DO ig=1,ngrid
178      unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
179                                                      ENDDO
180      DO k=2,nlay
181                                                      DO ig=1,ngrid
182        unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
183                                                     ENDDO
184      ENDDO
185                                                      DO ig=1,ngrid
186      unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
187                                                     ENDDO
188c
189c.......................................................................
190
191      do k=2,klev
192                                                          do ig=1,ngrid
193         dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
194         m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)
195     s             /(dz(ig,k)*dz(ig,k))
196         dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
197         n2(ig,k)=g*2.*dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
198c        n2(ig,k)=0.
199         ri=n2(ig,k)/max(m2(ig,k),1.e-10)
200         if (ri.lt.ric) then
201            rif(ig,k)=frif(ri)
202         else
203            rif(ig,k)=rifc
204         endif
205         if(rif(ig,k).lt.0.16) then
206            alpha(ig,k)=falpha(rif(ig,k))
207            sm(ig,k)=fsm(rif(ig,k))
208         else
209            alpha(ig,k)=1.12
210            sm(ig,k)=0.085
211         endif
212         zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
213c     print*,'RIF L=',k,rif(ig,k),ri*alpha(ig,k)
214
215
216                                                          enddo
217      enddo
218
219
220c====================================================================
221c   Au premier appel, on determine l et q2 de facon iterative.
222c iterration pour determiner la longueur de melange
223
224
225      if (first.or.iflag_pbl.eq.6) then
226                                                          do ig=1,ngrid
227      l0(ig)=10.
228                                                          enddo
229      do k=2,klev-1
230                                                          do ig=1,ngrid
231        l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
232                                                          enddo
233      enddo
234
235      do iter=1,10
236                                                          do ig=1,ngrid
237         sq(ig)=1.e-10
238         sqz(ig)=1.e-10
239                                                          enddo
240         do k=2,klev-1
241                                                          do ig=1,ngrid
242           q2(ig,k)=l(ig,k)**2*zz(ig,k)
243           l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
244           zq=sqrt(q2(ig,k))
245           sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
246           sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
247                                                          enddo
248         enddo
249                                                          do ig=1,ngrid
250         l0(ig)=0.2*sqz(ig)/sq(ig)
251c        l0(ig)=30.
252                                                          enddo
253c      print*,'ITER=',iter,'  L0=',l0
254
255      enddo
256
257c     print*,'Fin de l initialisation de q2 et l0'
258
259      endif ! first
260
261c====================================================================
262c  Calcul de la longueur de melange.
263c====================================================================
264
265c   Mise a jour de l0
266                                                          do ig=1,ngrid
267      sq(ig)=1.e-10
268      sqz(ig)=1.e-10
269                                                          enddo
270      do k=2,klev-1
271                                                          do ig=1,ngrid
272        zq=sqrt(q2(ig,k))
273        sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
274        sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
275                                                          enddo
276      enddo
277                                                          do ig=1,ngrid
278      l0(ig)=0.2*sqz(ig)/sq(ig)
279c        l0(ig)=30.
280                                                          enddo
281c      print*,'ITER=',iter,'  L0=',l0
282c   calcul de l(z)
283      do k=2,klev
284                                                          do ig=1,ngrid
285         l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
286         if(first) then
287           q2(ig,k)=l(ig,k)**2*zz(ig,k)
288         endif
289                                                          enddo
290      enddo
291
292c====================================================================
293c   Yamada 2.0
294c====================================================================
295      if (iflag_pbl.eq.6) then
296
297      do k=2,klev
298                                                          do ig=1,ngrid
299         q2(ig,k)=l(ig,k)**2*zz(ig,k)
300                                                          enddo
301      enddo
302
303
304      else if (iflag_pbl.eq.7) then
305c====================================================================
306c   Yamada 2.Fournier
307c====================================================================
308
309c  Calcul de l,  km, au pas precedent
310      do k=2,klev
311                                                          do ig=1,ngrid
312c        print*,'SMML=',sm(ig,k),l(ig,k)
313         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
314         kmpre(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
315         mpre(ig,k)=sqrt(m2(ig,k))
316c        print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
317                                                          enddo
318      enddo
319
320      do k=2,klev-1
321                                                          do ig=1,ngrid
322        m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
323        mcstat=sqrt(m2cstat)
324
325c        print*,'M2 L=',k,mpre(ig,k),mcstat
326c
327c  -----{puis on ecrit la valeur de q qui annule l'equation de m
328c        supposee en q3}
329c
330        IF (k.eq.2) THEN
331          kmcstat=1.E+0 / mcstat
332     &    *( unsdz(ig,k)*kmpre(ig,k+1)
333     &                        *mpre(ig,k+1)
334     &      +unsdz(ig,k-1)
335     &              *cd(ig)
336     &              *( sqrt(u(ig,3)**2+v(ig,3)**2)
337     &                -mcstat/unsdzdec(ig,k)
338     &                -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
339     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
340        ELSE
341          kmcstat=1.E+0 / mcstat
342     &    *( unsdz(ig,k)*kmpre(ig,k+1)
343     &                        *mpre(ig,k+1)
344     &      +unsdz(ig,k-1)*kmpre(ig,k-1)
345     &                          *mpre(ig,k-1) )
346     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
347        ENDIF
348c       print*,'T2 L=',k,tmp2
349        tmp2=kmcstat
350     &      /( sm(ig,k)/q2(ig,k) )
351     &      /l(ig,k)
352        q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
353c       print*,'Q2 L=',k,q2(ig,k)
354c
355                                                          enddo
356      enddo
357
358      else if (iflag_pbl.ge.8) then
359c====================================================================
360c   Yamada 2.5 a la Didi
361c====================================================================
362
363
364c  Calcul de l,  km, au pas precedent
365      do k=2,klev
366                                                          do ig=1,ngrid
367c        print*,'SMML=',sm(ig,k),l(ig,k)
368         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
369         if (delta(ig,k).lt.1.e-20) then
370c     print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
371            delta(ig,k)=1.e-20
372         endif
373         km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
374         aa0=
375     s   (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
376         aa1=
377     s   (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
378c abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
379         aa(ig,k)=aa1*dt/(delta(ig,k)*l(ig,k))
380c     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
381         qpre=sqrt(q2(ig,k))
382         if (iflag_pbl.eq.8 ) then
383            if (aa(ig,k).gt.0.) then
384               q2(ig,k)=(qpre+aa(ig,k)*qpre*qpre)**2
385            else
386               q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
387            endif
388         else ! iflag_pbl=9
389            if (aa(ig,k)*qpre.gt.0.9) then
390               q2(ig,k)=(qpre*10.)**2
391            else
392               q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
393            endif
394         endif
395         q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
396c     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
397                                                          enddo
398      enddo
399
400      endif ! Fin du cas 8
401
402c     print*,'OK8'
403
404c====================================================================
405c   Calcul des coefficients de m�ange
406c====================================================================
407      do k=2,klev
408c     print*,'k=',k
409                                                          do ig=1,ngrid
410cabde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
411         zq=sqrt(q2(ig,k))
412         km(ig,k)=l(ig,k)*zq*sm(ig,k)
413         kn(ig,k)=km(ig,k)*alpha(ig,k)
414         kq(ig,k)=l(ig,k)*zq*0.2
415c     print*,'KML=',km(ig,k),kn(ig,k)
416                                                          enddo
417      enddo
418
419c     if (iflag_pbl.ge.7..and.0.eq.1) then
420c        q2(:,1)=q2(:,2)
421c        call vdif_q2(dt,g,rconst,plev,temp,kq,q2)
422c     endif
423
424c   Traitement des cas noctrunes avec l'introduction d'une longueur
425c   minilale.
426
427c====================================================================
428c   Traitement particulier pour les cas tres stables.
429c   D'apres Holtslag Boville.
430
431      if (prt_level>1) THEN
432       print*,'YAMADA4 0'
433      endif !(prt_level>1) THEN
434                                                          do ig=1,ngrid
435      coriol(ig)=1.e-4
436      pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
437                                                          enddo
438
439!      print*,'pblhmin ',pblhmin
440CTest a remettre 21 11 02
441c test abd 13 05 02      if(0.eq.1) then
442      if(1.eq.1) then
443      do k=2,klev
444         do ig=1,ngrid
445            if (teta(ig,2).gt.teta(ig,1)) then
446               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
447               kmin=kap*zlev(ig,k)*qmin
448            else
449               kmin=-1. ! kmin n'est utilise que pour les SL stables.
450            endif
451            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
452c               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
453c     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
454               kn(ig,k)=kmin
455               km(ig,k)=kmin
456               kq(ig,k)=kmin
457c   la longueur de melange est suposee etre l= kap z
458c   K=l q Sm d'ou q2=(K/l Sm)**2
459               q2(ig,k)=(qmin/sm(ig,k))**2
460            endif
461         enddo
462      enddo
463      endif
464
465      if (prt_level>1) THEN
466       print*,'YAMADA4 1'
467      endif !(prt_level>1) THEN
468c   Diagnostique pour stokage
469
470      if(1.eq.0)then
471      rino=rif
472      smyam(1:ngrid,1)=0.
473      styam(1:ngrid,1)=0.
474      lyam(1:ngrid,1)=0.
475      knyam(1:ngrid,1)=0.
476      w2yam(1:ngrid,1)=0.
477      t2yam(1:ngrid,1)=0.
478
479      smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)
480      styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev)
481      lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev)
482      knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev)
483
484c   Estimations de w'2 et T'2 d'apres Abdela et McFarlane
485
486      w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24
487     s    +lyam(1:ngrid,2:klev)*5.17*kn(1:ngrid,2:klev)
488     s    *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev))
489
490      t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev)
491     s    *dtetadz(1:ngrid,2:klev)**2
492     s    /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev)
493      endif
494
495c     print*,'OKFIN'
496      first=.false.
497      return
498      end
Note: See TracBrowser for help on using the repository browser.