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

Diff of /trunk/phylmd/yamada4.f

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

revision 117 by guez, Wed Jul 2 18:39:15 2014 UTC revision 118 by guez, Thu Dec 18 17:30:24 2014 UTC
# Line 2  module yamada4_m Line 2  module yamada4_m
2    
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    real, parameter, private:: kap = 0.4    private
6      public yamada4
7      real, parameter:: kap = 0.4
8    
9  contains  contains
10    
11    SUBROUTINE yamada4(ngrid, dt, g, zlev, zlay, u, v, teta, cd, q2, km, kn, &    SUBROUTINE yamada4(ngrid, dt, g, zlev, zlay, u, v, teta, cd, q2, km, kn, kq, &
12         kq, ustar, iflag_pbl)         ustar, iflag_pbl)
13    
14      ! From LMDZ4/libf/phylmd/yamada4.F, version 1.1 2004/06/22 11:45:36      ! From LMDZ4/libf/phylmd/yamada4.F, version 1.1 2004/06/22 11:45:36
15    
16      USE dimphy, ONLY : klev, klon      use nr_util, only: assert
17        USE dimphy, ONLY: klev
18    
19      integer ngrid      integer, intent(in):: ngrid
20      REAL, intent(in):: dt ! pas de temps      REAL, intent(in):: dt ! pas de temps
21      real, intent(in):: g      real, intent(in):: g
22    
23      REAL zlev(klon, klev+1)      REAL zlev(ngrid, klev+1)
24      ! altitude à chaque niveau (interface inférieure de la couche de      ! altitude à chaque niveau (interface inférieure de la couche de
25      ! même indice)      ! même indice)
26    
27      REAL zlay(klon, klev) ! altitude au centre de chaque couche      REAL zlay(ngrid, klev) ! altitude au centre de chaque couche
28    
29      REAL u(klon, klev), v(klon, klev)      REAL u(ngrid, klev), v(ngrid, klev)
30      ! vitesse au centre de chaque couche (en entrée : la valeur au      ! vitesse au centre de chaque couche (en entrée : la valeur au
31      ! début du pas de temps)      ! début du pas de temps)
32    
33      REAL teta(klon, klev)      REAL, intent(in): teta(ngrid, klev)
34      ! température potentielle au centre de chaque couche (en entrée :      ! température potentielle au centre de chaque couche (en entrée :
35      ! la valeur au début du pas de temps)      ! la valeur au début du pas de temps)
36    
37      REAL, intent(in):: cd(:) ! (ngrid) cdrag, valeur au début du pas de temps      REAL, intent(in):: cd(:) ! (ngrid) cdrag, valeur au début du pas de temps
38    
39      REAL, intent(inout):: q2(klon, klev+1)      REAL, intent(inout):: q2(ngrid, klev+1)
40      ! $q^2$ au bas de chaque couche      ! $q^2$ au bas de chaque couche
41      ! En entrée : la valeur au début du pas de temps ; en sortie : la      ! En entrée : la valeur au début du pas de temps ; en sortie : la
42      ! valeur à la fin du pas de temps.      ! valeur à la fin du pas de temps.
43    
44      REAL km(klon, klev+1)      REAL km(ngrid, klev+1)
45      ! diffusivité turbulente de quantité de mouvement (au bas de      ! diffusivité turbulente de quantité de mouvement (au bas de
46      ! chaque couche) (en sortie : la valeur à la fin du pas de temps)      ! chaque couche) (en sortie : la valeur à la fin du pas de temps)
47    
48      REAL kn(klon, klev+1)      REAL kn(ngrid, klev+1)
49      ! diffusivité turbulente des scalaires (au bas de chaque couche)      ! diffusivité turbulente des scalaires (au bas de chaque couche)
50      ! (en sortie : la valeur à la fin du pas de temps)      ! (en sortie : la valeur à la fin du pas de temps)
51    
52      REAL kq(klon, klev+1)      REAL kq(ngrid, klev+1)
53      real ustar(klon)      real ustar(ngrid)
54    
55      integer iflag_pbl      integer iflag_pbl
56      ! iflag_pbl doit valoir entre 6 et 9      ! iflag_pbl doit valoir entre 6 et 9
# Line 59  contains Line 62  contains
62    
63      ! Local:      ! Local:
64    
65      real kmin, qmin, pblhmin(klon), coriol(klon)      real kmin, qmin, pblhmin(ngrid), coriol(ngrid)
66      real qpre      real qpre
67      REAL unsdz(klon, klev)      REAL unsdz(ngrid, klev)
68      REAL unsdzdec(klon, klev+1)      REAL unsdzdec(ngrid, klev+1)
69      REAL kmpre(klon, klev+1), tmp2      REAL kmpre(ngrid, klev+1), tmp2
70      REAL mpre(klon, klev+1)      REAL mpre(ngrid, klev+1)
71      real delta(klon, klev+1)      real delta(ngrid, klev+1)
72      real aa(klon, klev+1), aa0, aa1      real aa(ngrid, klev+1), aa0, aa1
     integer, PARAMETER:: nlay = klev  
73      integer, PARAMETER:: nlev = klev+1      integer, PARAMETER:: nlev = klev+1
74      logical:: first = .true.      logical:: first = .true.
75      integer:: ipas = 0      integer:: ipas = 0
76      integer ig, k      integer ig, k
77      real ri      real ri
78      real rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)      real rif(ngrid, klev+1), sm(ngrid, klev+1), alpha(ngrid, klev)
79      real m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)      real m2(ngrid, klev+1), dz(ngrid, klev+1), zq, n2(ngrid, klev+1)
80      real dtetadz(klon, klev+1)      real dtetadz(ngrid, klev+1)
81      real m2cstat, mcstat, kmcstat      real m2cstat, mcstat, kmcstat
82      real l(klon, klev+1)      real l(ngrid, klev+1)
83      real, save:: l0(klon)      real, save:: l0(ngrid)
84      real sq(klon), sqz(klon), zz(klon, klev+1)      real sq(ngrid), sqz(ngrid), zz(ngrid, klev+1)
85      integer iter      integer iter
86      real:: ric = 0.195, rifc = 0.191, b1 = 16.6, kap = 0.4      real:: ric = 0.195, rifc = 0.191, b1 = 16.6, kap = 0.4
87      real rino(klon, klev+1), smyam(klon, klev), styam(klon, klev)      real rino(ngrid, klev+1), smyam(ngrid, klev), styam(ngrid, klev)
88      real lyam(klon, klev), knyam(klon, klev)      real lyam(ngrid, klev), knyam(ngrid, klev)
89    
90      !-----------------------------------------------------------------------      !-----------------------------------------------------------------------
91    
92      if (.not. (iflag_pbl >= 6 .and. iflag_pbl <= 9)) then      call assert(iflag_pbl >= 6 .and. iflag_pbl <= 9, "yamada4")
        print *, 'probleme de coherence dans appel a MY'  
        stop 1  
     endif  
93    
94      ipas = ipas+1      ipas = ipas+1
95    
96      ! les increments verticaux      ! les increments verticaux
97      DO ig = 1, ngrid      DO ig = 1, ngrid
98         ! alerte: zlev n'est pas declare a nlev         ! alerte: zlev n'est pas declare a nlev
99         zlev(ig, nlev) = zlay(ig, nlay) +(zlay(ig, nlay) - zlev(ig, nlev-1))         zlev(ig, nlev) = zlay(ig, klev) +(zlay(ig, klev) - zlev(ig, nlev-1))
100      ENDDO      ENDDO
101    
102      DO k = 1, nlay      DO k = 1, klev
103         DO ig = 1, ngrid         DO ig = 1, ngrid
104            unsdz(ig, k) = 1.E+0/(zlev(ig, k+1)-zlev(ig, k))            unsdz(ig, k) = 1.E+0/(zlev(ig, k+1)-zlev(ig, k))
105         ENDDO         ENDDO
106      ENDDO      ENDDO
107    
108      DO ig = 1, ngrid      DO ig = 1, ngrid
109         unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1))         unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1))
110      ENDDO      ENDDO
111      DO k = 2, nlay  
112        DO k = 2, klev
113         DO ig = 1, ngrid         DO ig = 1, ngrid
114            unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1))            unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1))
115         ENDDO         ENDDO
116      ENDDO      ENDDO
117    
118      DO ig = 1, ngrid      DO ig = 1, ngrid
119         unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig, nlay+1)-zlay(ig, nlay))         unsdzdec(ig, klev+1) = 1.E+0/(zlev(ig, klev+1)-zlay(ig, klev))
120      ENDDO      ENDDO
121    
122      do k = 2, klev      do k = 2, klev
# Line 309  contains Line 311  contains
311    
312      print *, 'pblhmin ', pblhmin      print *, 'pblhmin ', pblhmin
313      do k = 2, klev      do k = 2, klev
314         do ig = 1, klon         do ig = 1, ngrid
315            if (teta(ig, 2).gt.teta(ig, 1)) then            if (teta(ig, 2).gt.teta(ig, 1)) then
316               qmin = ustar(ig)*(max(1.-zlev(ig, k)/pblhmin(ig), 0.))**2               qmin = ustar(ig)*(max(1.-zlev(ig, k)/pblhmin(ig), 0.))**2
317               kmin = kap*zlev(ig, k)*qmin               kmin = kap*zlev(ig, k)*qmin
# Line 332  contains Line 334  contains
334      rino = rif      rino = rif
335      smyam(:, 1:klev) = sm(:, 1:klev)      smyam(:, 1:klev) = sm(:, 1:klev)
336      styam = sm(:, 1:klev)*alpha(:, 1:klev)      styam = sm(:, 1:klev)*alpha(:, 1:klev)
337      lyam(1:klon, 1:klev) = l(:, 1:klev)      lyam(1:ngrid, 1:klev) = l(:, 1:klev)
     knyam(1:klon, 1:klev) = kn(:, 1:klev)  
338    
339      first = .false.      first = .false.
340    
Legend:
Removed from v.117  
changed lines
  Added in v.118
  ViewVC Help
Powered by ViewVC 1.1.21