/[lmdze]/trunk/phylmd/Interface_surf/yamada4.f
ViewVC logotype

Diff of /trunk/phylmd/Interface_surf/yamada4.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/yamada4.f revision 12 by guez, Mon Jul 21 16:05:07 2008 UTC trunk/phylmd/Interface_surf/yamada4.f revision 303 by guez, Thu Sep 6 14:25:07 2018 UTC
# Line 1  Line 1 
1  !  module yamada4_m
2  ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/yamada4.F,v 1.1 2004/06/22 11:45:36 lmdzadmin Exp $  
3  !    IMPLICIT NONE
4        SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp  
5       s   ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar    private
6       s   ,iflag_pbl)    public yamada4
7        use dimens_m    real, parameter:: kap = 0.4
8        use dimphy  
9        IMPLICIT NONE  contains
10  c.......................................................................  
11  c.......................................................................    SUBROUTINE yamada4(zlev, zlay, u, v, teta, q2, coefm, coefh, ustar)
12  c  
13  c dt : pas de temps      ! From LMDZ4/libf/phylmd/yamada4.F, version 1.1 2004/06/22 11:45:36
14  c g  : g  
15  c zlev : altitude a chaque niveau (interface inferieure de la couche      ! Library:
16  c        de meme indice)      use nr_util, only: assert, assert_eq
17  c zlay : altitude au centre de chaque couche  
18  c u,v : vitesse au centre de chaque couche      use comconst, only: dtphys
19  c       (en entree : la valeur au debut du pas de temps)      USE conf_phys_m, ONLY: iflag_pbl
20  c teta : temperature potentielle au centre de chaque couche      USE dimphy, ONLY: klev
21  c        (en entree : la valeur au debut du pas de temps)      USE suphec_m, ONLY: rg
22  c cd : cdrag  
23  c      (en entree : la valeur au debut du pas de temps)      REAL, intent(in):: zlev(:, :) ! (knon, klev + 1)
24  c q2 : $q^2$ au bas de chaque couche      ! altitude \`a chaque niveau (interface inf\'erieure de la couche de
25  c      (en entree : la valeur au debut du pas de temps)      ! m\^eme indice)
26  c      (en sortie : la valeur a la fin du pas de temps)  
27  c km : diffusivite turbulente de quantite de mouvement (au bas de chaque      REAL, intent(in):: zlay(:, :) ! (knon, klev) altitude au centre de
28  c      couche)                                    ! chaque couche
29  c      (en sortie : la valeur a la fin du pas de temps)  
30  c kn : diffusivite turbulente des scalaires (au bas de chaque couche)      REAL, intent(in):: u(:, :), v(:, :) ! (knon, klev)
31  c      (en sortie : la valeur a la fin du pas de temps)      ! vitesse au centre de chaque couche (en entr\'ee : la valeur au
32  c      ! d\'ebut du pas de temps)
33  c  iflag_pbl doit valoir entre 6 et 9  
34  c      l=6, on prend  systematiquement une longueur d'equilibre      REAL, intent(in):: teta(:, :) ! (knon, klev)
35  c    iflag_pbl=6 : MY 2.0      ! temp\'erature potentielle au centre de chaque couche (en entr\'ee :
36  c    iflag_pbl=7 : MY 2.0.Fournier      ! la valeur au d\'ebut du pas de temps)
37  c    iflag_pbl=8 : MY 2.5  
38  c    iflag_pbl=9 : un test ?      REAL, intent(inout):: q2(:, :) ! (knon, klev + 1)
39        ! $q^2$ au bas de chaque couche
40  c.......................................................................      ! En entr\'ee : la valeur au d\'ebut du pas de temps ; en sortie : la
41        REAL, intent(in):: dt      ! valeur \`a la fin du pas de temps.
42        real g,rconst  
43        real plev(klon,klev+1),temp(klon,klev)      REAL, intent(out):: coefm(:, 2:) ! (knon, 2:klev)
44        real ustar(klon)      ! diffusivit\'e turbulente de quantit\'e de mouvement (au bas de
45        real kmin,qmin,pblhmin(klon),coriol(klon)      ! chaque couche) (en sortie : la valeur \`a la fin du pas de temps)
46        REAL zlev(klon,klev+1)  
47        REAL zlay(klon,klev)      REAL, intent(out):: coefh(:, 2:) ! (knon, 2:klev)
48        REAL u(klon,klev)      ! diffusivit\'e turbulente des scalaires (au bas de chaque couche)
49        REAL v(klon,klev)      ! (en sortie : la valeur \`a la fin du pas de temps)
50        REAL teta(klon,klev)  
51        REAL cd(klon)      real, intent(in):: ustar(:) ! (knon)
52        REAL q2(klon,klev+1),qpre  
53        REAL unsdz(klon,klev)      ! Local:
54        REAL unsdzdec(klon,klev+1)      
55        integer knon
56        REAL km(klon,klev+1)      real kmin, qmin
57        REAL kmpre(klon,klev+1),tmp2      real pblhmin(size(ustar)), coriol(size(ustar)) ! (knon)
58        REAL mpre(klon,klev+1)      real qpre
59        REAL kn(klon,klev+1)      REAL unsdz(size(zlay, 1), size(zlay, 2)) ! (knon, klev)
60        REAL kq(klon,klev+1)      REAL unsdzdec(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
61        real ff(klon,klev+1),delta(klon,klev+1)      real delta(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
62        real aa(klon,klev+1),aa0,aa1      real aa(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
63        integer iflag_pbl,ngrid      logical:: first = .true.
64        integer ig, k
65        real ri
66        integer nlay,nlev      real, dimension(size(zlev, 1), size(zlev, 2)):: rif, sm ! (knon, klev + 1)
67        PARAMETER (nlay=klev)      real alpha(size(zlay, 1), size(zlay, 2)) ! (knon, klev)
68        PARAMETER (nlev=klev+1)  
69        real, dimension(size(zlev, 1), size(zlev, 2)):: m2, dz, n2
70        logical first      ! (knon, klev + 1)
71        integer ipas      
72        save first,ipas      real zq
73        data first,ipas/.true.,0/      real dtetadz(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
74        real l(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
75        real l0(size(ustar)) ! (knon)
76        integer ig,k      real sq(size(ustar)), sqz(size(ustar)) ! (knon)
77        real zz(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
78        integer iter
79        real ri,zrif,zalpha,zsm,zsn      real:: ric = 0.195, rifc = 0.191, b1 = 16.6
80        real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)  
81        !-----------------------------------------------------------------------
82        real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)  
83        real dtetadz(klon,klev+1)      call assert(any(iflag_pbl == [6, 8, 9]), "yamada4 iflag_pbl")
84        real m2cstat,mcstat,kmcstat      knon = assert_eq([size(zlev, 1), size(zlay, 1), size(u, 1), size(v, 1), &
85        real l(klon,klev+1),l0(klon)           size(teta, 1), size(ustar), size(q2, 1), size(coefm, 1), &
86        save l0           size(coefh, 1)], "yamada4 knon")
87        call assert(klev == [size(zlev, 2) - 1, size(zlay, 2), size(u, 2), &
88        real sq(klon),sqz(klon),zz(klon,klev+1)           size(v, 2), size(teta, 2), size(q2, 2) - 1, size(coefm, 2) + 1, &
89        integer iter           size(coefh, 2) + 1], "yamada4 klev")
90    
91        real ric,rifc,b1,kap      ! les increments verticaux
92        save ric,rifc,b1,kap  
93        data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/      DO k = 1, klev
94           DO ig = 1, knon
95        real frif,falpha,fsm            unsdz(ig, k) = 1.E+0/(zlev(ig, k + 1)-zlev(ig, k))
96        real fl,zzz,zl0,zq2,zn2         ENDDO
97        ENDDO
98        real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)  
99       s  ,lyam(klon,klev),knyam(klon,klev)      DO ig = 1, knon
100       s  ,w2yam(klon,klev),t2yam(klon,klev)         unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1))
101        common/pbldiag/rino,smyam,styam,lyam,knyam,w2yam,t2yam      ENDDO
102    
103        frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))      DO k = 2, klev
104        falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)         DO ig = 1, knon
105        fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))            unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1))
106        fl(zzz,zl0,zq2,zn2)=         ENDDO
107       s     max(min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))      ENDDO
108       s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)  
109        DO ig = 1, knon
110        if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then         unsdzdec(ig, klev + 1) = 1.E+0/(zlev(ig, klev + 1)-zlay(ig, klev))
111             stop'probleme de coherence dans appel a MY'      ENDDO
112        endif  
113        do k = 2, klev
114        ipas=ipas+1         do ig = 1, knon
115        if (0.eq.1.and.first) then            dz(ig, k) = zlay(ig, k)-zlay(ig, k-1)
116        do ig=1,1000            m2(ig, k) = ((u(ig, k)-u(ig, k-1))**2 + (v(ig, k)-v(ig, k-1))**2) &
117           ri=(ig-800.)/500.                 /(dz(ig, k)*dz(ig, k))
118           if (ri.lt.ric) then            dtetadz(ig, k) = (teta(ig, k)-teta(ig, k-1))/dz(ig, k)
119              zrif=frif(ri)            n2(ig, k) = rg*2.*dtetadz(ig, k)/(teta(ig, k-1) + teta(ig, k))
120           else            ri = n2(ig, k)/max(m2(ig, k), 1.e-10)
121              zrif=rifc            if (ri.lt.ric) then
122           endif               rif(ig, k) = frif(ri)
123           if(zrif.lt.0.16) then            else
124              zalpha=falpha(zrif)               rif(ig, k) = rifc
125              zsm=fsm(zrif)            endif
126           else            if (rif(ig, k).lt.0.16) then
127              zalpha=1.12               alpha(ig, k) = falpha(rif(ig, k))
128              zsm=0.085               sm(ig, k) = fsm(rif(ig, k))
129           endif            else
130  c     print*,ri,rif,zalpha,zsm               alpha(ig, k) = 1.12
131        enddo               sm(ig, k) = 0.085
132        endif            endif
133              zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig, k))*sm(ig, k)
134  c.......................................................................         enddo
135  c  les increments verticaux      enddo
136  c.......................................................................  
137  c      ! Au premier appel, on d\'etermine l et q2 de fa\ccon it\'erative.
138  c!!!!! allerte !!!!!c      ! It\'eration pour d\'eterminer la longueur de m\'elange
139  c!!!!! zlev n'est pas declare a nlev !!!!!c  
140  c!!!!! ---->      if (first .or. iflag_pbl == 6) then
141                                                        DO ig=1,ngrid         do ig = 1, knon
142              zlev(ig,nlev)=zlay(ig,nlay)            l0(ig) = 10.
143       &             +( zlay(ig,nlay) - zlev(ig,nlev-1) )         enddo
144                                                        ENDDO         do k = 2, klev-1
145  c!!!!! <----            do ig = 1, knon
146  c!!!!! allerte !!!!!c               l(ig, k) = l0(ig) * kap * zlev(ig, k) &
147  c                    / (kap * zlev(ig, k) + l0(ig))
148        DO k=1,nlay            enddo
149                                                        DO ig=1,ngrid         enddo
150          unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))  
151                                                        ENDDO         do iter = 1, 10
152        ENDDO            do ig = 1, knon
153                                                        DO ig=1,ngrid               sq(ig) = 1e-10
154        unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))               sqz(ig) = 1e-10
155                                                        ENDDO            enddo
156        DO k=2,nlay            do k = 2, klev-1
157                                                        DO ig=1,ngrid               do ig = 1, knon
158          unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))                  q2(ig, k) = l(ig, k)**2 * zz(ig, k)
159                                                       ENDDO                  l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k))
160        ENDDO                  zq = sqrt(q2(ig, k))
161                                                        DO ig=1,ngrid                  sqz(ig) = sqz(ig) + zq * zlev(ig, k) &
162        unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))                       * (zlay(ig, k) - zlay(ig, k-1))
163                                                       ENDDO                  sq(ig) = sq(ig) + zq * (zlay(ig, k) - zlay(ig, k-1))
164  c               enddo
165  c.......................................................................            enddo
166              do ig = 1, knon
167        do k=2,klev               l0(ig) = 0.2 * sqz(ig) / sq(ig)
168                                                            do ig=1,ngrid            enddo
169           dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)         enddo
170           m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)      endif
171       s             /(dz(ig,k)*dz(ig,k))  
172           dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)      ! Calcul de la longueur de melange.
173           n2(ig,k)=g*2.*dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))  
174  c        n2(ig,k)=0.      ! Mise a jour de l0
175           ri=n2(ig,k)/max(m2(ig,k),1.e-10)      do ig = 1, knon
176           if (ri.lt.ric) then         sq(ig) = 1.e-10
177              rif(ig,k)=frif(ri)         sqz(ig) = 1.e-10
178           else      enddo
179              rif(ig,k)=rifc      do k = 2, klev-1
180           endif         do ig = 1, knon
181           if(rif(ig,k).lt.0.16) then            zq = sqrt(q2(ig, k))
182              alpha(ig,k)=falpha(rif(ig,k))            sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig, k)-zlay(ig, k-1))
183              sm(ig,k)=fsm(rif(ig,k))            sq(ig) = sq(ig) + zq*(zlay(ig, k)-zlay(ig, k-1))
184           else         enddo
185              alpha(ig,k)=1.12      enddo
186              sm(ig,k)=0.085      do ig = 1, knon
187           endif         l0(ig) = 0.2*sqz(ig)/sq(ig)
188           zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)      enddo
189  c     print*,'RIF L=',k,rif(ig,k),ri*alpha(ig,k)      ! calcul de l(z)
190        do k = 2, klev
191           do ig = 1, knon
192                                                            enddo            l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k))
193        enddo            if (first) then
194                 q2(ig, k) = l(ig, k)**2 * zz(ig, k)
195              endif
196  c====================================================================         enddo
197  c   Au premier appel, on determine l et q2 de facon iterative.      enddo
198  c iterration pour determiner la longueur de melange  
199        if (iflag_pbl == 6) then
200           ! Yamada 2.0
201        if (first.or.iflag_pbl.eq.6) then         do k = 2, klev
202                                                            do ig=1,ngrid            do ig = 1, knon
203        l0(ig)=10.               q2(ig, k) = l(ig, k)**2 * zz(ig, k)
204                                                            enddo            enddo
205        do k=2,klev-1         enddo
206                                                            do ig=1,ngrid      else if (iflag_pbl >= 8) then
207          l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))         ! Yamada 2.5 a la Didi
208                                                            enddo  
209        enddo         ! Calcul de l, coefm, au pas precedent
210           do k = 2, klev
211        do iter=1,10            do ig = 1, knon
212                                                            do ig=1,ngrid               delta(ig, k) = q2(ig, k)/(l(ig, k)**2*sm(ig, k))
213           sq(ig)=1.e-10               if (delta(ig, k).lt.1.e-20) then
214           sqz(ig)=1.e-10                  delta(ig, k) = 1.e-20
215                                                            enddo               endif
216           do k=2,klev-1               coefm(ig, k) = l(ig, k)*sqrt(q2(ig, k))*sm(ig, k)
217                                                            do ig=1,ngrid               aa(ig, k) = (m2(ig, k)*(1.-rif(ig, k))-delta(ig, k)/b1)*dtphys/(delta(ig, k)*l(ig, k))
218             q2(ig,k)=l(ig,k)**2*zz(ig,k)               qpre = sqrt(q2(ig, k))
219             l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))               if (iflag_pbl == 8) then
220             zq=sqrt(q2(ig,k))                  if (aa(ig, k).gt.0.) then
221             sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))                     q2(ig, k) = (qpre + aa(ig, k)*qpre*qpre)**2
222             sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))                  else
223                                                            enddo                     q2(ig, k) = (qpre/(1.-aa(ig, k)*qpre))**2
224           enddo                  endif
225                                                            do ig=1,ngrid               else
226           l0(ig)=0.2*sqz(ig)/sq(ig)                  ! iflag_pbl = 9
227  c        l0(ig)=30.                  if (aa(ig, k)*qpre.gt.0.9) then
228                                                            enddo                     q2(ig, k) = (qpre*10.)**2
229  c      print*,'ITER=',iter,'  L0=',l0                  else
230                       q2(ig, k) = (qpre/(1.-aa(ig, k)*qpre))**2
231        enddo                  endif
232                 endif
233  c     print*,'Fin de l initialisation de q2 et l0'               q2(ig, k) = min(max(q2(ig, k), 1.e-10), 1.e4)
234              enddo
235        endif ! first         enddo
236        endif
237  c====================================================================  
238  c  Calcul de la longueur de melange.      ! Calcul des coefficients de m\'elange
239  c====================================================================      do k = 2, klev
240           do ig = 1, knon
241  c   Mise a jour de l0            zq = sqrt(q2(ig, k))
242                                                            do ig=1,ngrid            coefm(ig, k) = l(ig, k)*zq*sm(ig, k)
243        sq(ig)=1.e-10            coefh(ig, k) = coefm(ig, k)*alpha(ig, k)
244        sqz(ig)=1.e-10         enddo
245                                                            enddo      enddo
246        do k=2,klev-1  
247                                                            do ig=1,ngrid      ! Traitement des cas noctrunes avec l'introduction d'une longueur
248          zq=sqrt(q2(ig,k))      ! minilale.
249          sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))  
250          sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))      ! Traitement particulier pour les cas tres stables.
251                                                            enddo      ! D'apres Holtslag Boville.
252        enddo  
253                                                            do ig=1,ngrid      do ig = 1, knon
254        l0(ig)=0.2*sqz(ig)/sq(ig)         coriol(ig) = 1.e-4
255  c        l0(ig)=30.         pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546e-5)
256                                                            enddo      enddo
257  c      print*,'ITER=',iter,'  L0=',l0  
258  c   calcul de l(z)      do k = 2, klev
259        do k=2,klev         do ig = 1, knon
260                                                            do ig=1,ngrid            if (teta(ig, 2).gt.teta(ig, 1)) then
261           l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))               qmin = ustar(ig)*(max(1.-zlev(ig, k)/pblhmin(ig), 0.))**2
262           if(first) then               kmin = kap*zlev(ig, k)*qmin
263             q2(ig,k)=l(ig,k)**2*zz(ig,k)            else
264           endif               kmin = -1. ! kmin n'est utilise que pour les SL stables.
265                                                            enddo            endif
266        enddo            if (coefh(ig, k).lt.kmin.or.coefm(ig, k).lt.kmin) then
267                 coefh(ig, k) = kmin
268  c====================================================================               coefm(ig, k) = kmin
269  c   Yamada 2.0               ! la longueur de melange est suposee etre l = kap z
270  c====================================================================               ! K = l q Sm d'ou q2 = (K/l Sm)**2
271        if (iflag_pbl.eq.6) then               q2(ig, k) = (qmin/sm(ig, k))**2
272              endif
273        do k=2,klev         enddo
274                                                            do ig=1,ngrid      enddo
275           q2(ig,k)=l(ig,k)**2*zz(ig,k)  
276                                                            enddo      first = .false.
277        enddo  
278      end SUBROUTINE yamada4
279    
280        else if (iflag_pbl.eq.7) then    !*******************************************************************
281  c====================================================================  
282  c   Yamada 2.Fournier    pure real function frif(ri)
283  c====================================================================  
284        real, intent(in):: ri
285  c  Calcul de l,  km, au pas precedent  
286        do k=2,klev      frif = 0.6588*(ri + 0.1776-sqrt(ri*ri-0.3221*ri + 0.03156))
287                                                            do ig=1,ngrid  
288  c        print*,'SMML=',sm(ig,k),l(ig,k)    end function frif
289           delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))  
290           kmpre(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)    !*******************************************************************
291           mpre(ig,k)=sqrt(m2(ig,k))  
292  c        print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)    pure real function falpha(ri)
293                                                            enddo  
294        enddo      real, intent(in):: ri
295    
296        do k=2,klev-1      falpha = 1.318*(0.2231-ri)/(0.2341-ri)
297                                                            do ig=1,ngrid  
298          m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)    end function falpha
299          mcstat=sqrt(m2cstat)  
300      !*******************************************************************
301  c        print*,'M2 L=',k,mpre(ig,k),mcstat  
302  c    pure real function fsm(ri)
303  c  -----{puis on ecrit la valeur de q qui annule l'equation de m  
304  c        supposee en q3}      real, intent(in):: ri
305  c  
306          IF (k.eq.2) THEN      fsm = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
307            kmcstat=1.E+0 / mcstat  
308       &    *( unsdz(ig,k)*kmpre(ig,k+1)    end function fsm
309       &                        *mpre(ig,k+1)  
310       &      +unsdz(ig,k-1)    !*******************************************************************
311       &              *cd(ig)  
312       &              *( sqrt(u(ig,3)**2+v(ig,3)**2)    pure real function fl(zzz, zl0, zq2, zn2)
313       &                -mcstat/unsdzdec(ig,k)  
314       &                -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)      real, intent(in):: zzz, zl0, zq2, zn2
315       &      /( unsdz(ig,k)+unsdz(ig,k-1) )  
316          ELSE      fl = max(min(zl0 * kap * zzz / (kap * zzz + zl0), &
317            kmcstat=1.E+0 / mcstat           0.5 * sqrt(zq2) / sqrt(max(zn2, 1e-10))), 1.)
318       &    *( unsdz(ig,k)*kmpre(ig,k+1)  
319       &                        *mpre(ig,k+1)    end function fl
320       &      +unsdz(ig,k-1)*kmpre(ig,k-1)  
321       &                          *mpre(ig,k-1) )  end module yamada4_m
      &      /( unsdz(ig,k)+unsdz(ig,k-1) )  
         ENDIF  
 c       print*,'T2 L=',k,tmp2  
         tmp2=kmcstat  
      &      /( sm(ig,k)/q2(ig,k) )  
      &      /l(ig,k)  
         q2(ig,k)=max(tmp2,1.e-12)**(2./3.)  
 c       print*,'Q2 L=',k,q2(ig,k)  
 c  
                                                           enddo  
       enddo  
   
       else if (iflag_pbl.ge.8) then  
 c====================================================================  
 c   Yamada 2.5 a la Didi  
 c====================================================================  
   
   
 c  Calcul de l,  km, au pas precedent  
       do k=2,klev  
                                                           do ig=1,ngrid  
 c        print*,'SMML=',sm(ig,k),l(ig,k)  
          delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))  
          if (delta(ig,k).lt.1.e-20) then  
 c     print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)  
             delta(ig,k)=1.e-20  
          endif  
          km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)  
          aa0=  
      s   (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)  
          aa1=  
      s   (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)  
 c abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)  
          aa(ig,k)=aa1*dt/(delta(ig,k)*l(ig,k))  
 c     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)  
          qpre=sqrt(q2(ig,k))  
          if (iflag_pbl.eq.8 ) then  
             if (aa(ig,k).gt.0.) then  
                q2(ig,k)=(qpre+aa(ig,k)*qpre*qpre)**2  
             else  
                q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2  
             endif  
          else ! iflag_pbl=9  
             if (aa(ig,k)*qpre.gt.0.9) then  
                q2(ig,k)=(qpre*10.)**2  
             else  
                q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2  
             endif  
          endif  
          q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)  
 c     print*,'Q2 L=',k,q2(ig,k),qpre*qpre  
                                                           enddo  
       enddo  
   
       endif ! Fin du cas 8  
   
 c     print*,'OK8'  
   
 c====================================================================  
 c   Calcul des coefficients de mélange  
 c====================================================================  
       do k=2,klev  
 c     print*,'k=',k  
                                                           do ig=1,ngrid  
 cabde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)  
          zq=sqrt(q2(ig,k))  
          km(ig,k)=l(ig,k)*zq*sm(ig,k)  
          kn(ig,k)=km(ig,k)*alpha(ig,k)  
          kq(ig,k)=l(ig,k)*zq*0.2  
 c     print*,'KML=',km(ig,k),kn(ig,k)  
                                                           enddo  
       enddo  
   
 c   Traitement des cas noctrunes avec l'introduction d'une longueur  
 c   minilale.  
   
 c====================================================================  
 c   Traitement particulier pour les cas tres stables.  
 c   D'apres Holtslag Boville.  
   
       print*,'YAMADA4 0'  
   
                                                           do ig=1,ngrid  
       coriol(ig)=1.e-4  
       pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)  
                                                           enddo  
   
        print*,'pblhmin ',pblhmin  
 CTest a remettre 21 11 02  
 c test abd 13 05 02      if(0.eq.1) then  
       if(1.eq.1) then  
       do k=2,klev  
          do ig=1,klon  
             if (teta(ig,2).gt.teta(ig,1)) then  
                qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2  
                kmin=kap*zlev(ig,k)*qmin  
             else  
                kmin=-1. ! kmin n'est utilise que pour les SL stables.  
             endif  
             if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then  
 c               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)  
 c     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)  
                kn(ig,k)=kmin  
                km(ig,k)=kmin  
                kq(ig,k)=kmin  
 c   la longueur de melange est suposee etre l= kap z  
 c   K=l q Sm d'ou q2=(K/l Sm)**2  
                q2(ig,k)=(qmin/sm(ig,k))**2  
             endif  
          enddo  
       enddo  
       endif  
   
       print*,'YAMADA4 1'  
 c   Diagnostique pour stokage  
   
       rino=rif  
       smyam(:,1:klev)=sm(:,1:klev)  
       styam=sm(:,1:klev)*alpha(:,1:klev)  
       lyam(1:klon,1:klev)=l(:,1:klev)  
       knyam(1:klon,1:klev)=kn(:,1:klev)  
   
 c   Estimations de w'2 et T'2 d'apres Abdela et McFarlane  
   
         if(1.eq.0)then  
       w2yam=q2(:,1:klev)*0.24  
      s    +lyam(:,1:klev)*5.17*kn(:,1:klev)*n2(:,1:klev)  
      s   /sqrt(q2(:,1:klev))  
   
       t2yam=9.1*kn(:,1:klev)*dtetadz(:,1:klev)**2/sqrt(q2(:,1:klev))  
      s  *lyam(:,1:klev)  
         endif  
   
 c     print*,'OKFIN'  
       first=.false.  
       return  
       end  

Legend:
Removed from v.12  
changed lines
  Added in v.303

  ViewVC Help
Powered by ViewVC 1.1.21