--- trunk/libf/phylmd/yamada4.f90 2011/07/01 15:00:48 47 +++ trunk/Sources/phylmd/yamada4.f 2016/03/11 18:47:26 178 @@ -2,58 +2,59 @@ IMPLICIT NONE - real, parameter:: kap = 0.4 private public yamada4 + real, parameter:: kap = 0.4 contains - SUBROUTINE yamada4(ngrid, dt, g, zlev, zlay, u, v, teta, cd, q2, km, kn, & - kq, ustar, iflag_pbl) + SUBROUTINE yamada4(ngrid, dt, g, zlev, zlay, u, v, teta, cd, q2, km, kn, kq, & + ustar, iflag_pbl) ! From LMDZ4/libf/phylmd/yamada4.F, version 1.1 2004/06/22 11:45:36 - USE dimphy, ONLY : klev, klon + use nr_util, only: assert + USE dimphy, ONLY: klev - integer ngrid + integer, intent(in):: ngrid REAL, intent(in):: dt ! pas de temps real, intent(in):: g - REAL zlev(klon, klev+1) - ! altitude à chaque niveau (interface inférieure de la couche de - ! même indice) + REAL zlev(ngrid, klev+1) + ! altitude \`a chaque niveau (interface inf\'erieure de la couche de + ! m\^eme indice) - REAL zlay(klon, klev) ! altitude au centre de chaque couche + REAL zlay(ngrid, klev) ! altitude au centre de chaque couche - REAL u(klon, klev), v(klon, klev) - ! vitesse au centre de chaque couche (en entrée : la valeur au - ! début du pas de temps) + REAL u(ngrid, klev), v(ngrid, klev) + ! vitesse au centre de chaque couche (en entr\'ee : la valeur au + ! d\'ebut du pas de temps) - REAL teta(klon, klev) - ! température potentielle au centre de chaque couche (en entrée : - ! la valeur au début du pas de temps) + REAL, intent(in):: teta(ngrid, klev) + ! temp\'erature potentielle au centre de chaque couche (en entr\'ee : + ! la valeur au d\'ebut du pas de temps) - REAL cd(klon) ! cdrag (en entrée : la valeur au début du pas de temps) + REAL, intent(in):: cd(:) ! (ngrid) cdrag, valeur au d\'ebut du pas de temps - REAL, intent(inout):: q2(klon, klev+1) + REAL, intent(inout):: q2(ngrid, klev+1) ! $q^2$ au bas de chaque couche - ! En entrée : la valeur au début du pas de temps ; en sortie : la - ! valeur à la fin du pas de temps. + ! En entr\'ee : la valeur au d\'ebut du pas de temps ; en sortie : la + ! valeur \`a la fin du pas de temps. - REAL km(klon, klev+1) - ! diffusivité turbulente de quantité de mouvement (au bas de - ! chaque couche) (en sortie : la valeur à la fin du pas de temps) + REAL km(ngrid, klev+1) + ! diffusivit\'e turbulente de quantit\'e de mouvement (au bas de + ! chaque couche) (en sortie : la valeur \`a la fin du pas de temps) - REAL kn(klon, klev+1) - ! diffusivité turbulente des scalaires (au bas de chaque couche) - ! (en sortie : la valeur à la fin du pas de temps) + REAL kn(ngrid, klev+1) + ! diffusivit\'e turbulente des scalaires (au bas de chaque couche) + ! (en sortie : la valeur \`a la fin du pas de temps) - REAL kq(klon, klev+1) - real ustar(klon) + REAL kq(ngrid, klev+1) + real ustar(ngrid) - integer iflag_pbl + integer, intent(in):: iflag_pbl ! iflag_pbl doit valoir entre 6 et 9 - ! l = 6, on prend systématiquement une longueur d'équilibre + ! l = 6, on prend syst\'ematiquement une longueur d'\'equilibre ! iflag_pbl = 6 : MY 2.0 ! iflag_pbl = 7 : MY 2.0.Fournier ! iflag_pbl = 8 : MY 2.5 @@ -61,62 +62,59 @@ ! Local: - real kmin, qmin, pblhmin(klon), coriol(klon) + real kmin, qmin, pblhmin(ngrid), coriol(ngrid) real qpre - REAL unsdz(klon, klev) - REAL unsdzdec(klon, klev+1) - REAL kmpre(klon, klev+1), tmp2 - REAL mpre(klon, klev+1) - real delta(klon, klev+1) - real aa(klon, klev+1), aa0, aa1 - integer, PARAMETER:: nlay = klev + REAL unsdz(ngrid, klev) + REAL unsdzdec(ngrid, klev+1) + REAL kmpre(ngrid, klev+1), tmp2 + REAL mpre(ngrid, klev+1) + real delta(ngrid, klev+1) + real aa(ngrid, klev+1), aa1 integer, PARAMETER:: nlev = klev+1 logical:: first = .true. integer:: ipas = 0 integer ig, k real ri - real rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev) - real m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1) - real dtetadz(klon, klev+1) + real rif(ngrid, klev+1), sm(ngrid, klev+1), alpha(ngrid, klev) + real m2(ngrid, klev+1), dz(ngrid, klev+1), zq, n2(ngrid, klev+1) + real dtetadz(ngrid, klev+1) real m2cstat, mcstat, kmcstat - real l(klon, klev+1) - real, save:: l0(klon) - real sq(klon), sqz(klon), zz(klon, klev+1) + real l(ngrid, klev+1) + real l0(ngrid) + real sq(ngrid), sqz(ngrid), zz(ngrid, klev+1) integer iter real:: ric = 0.195, rifc = 0.191, b1 = 16.6, kap = 0.4 - real rino(klon, klev+1), smyam(klon, klev), styam(klon, klev) - real lyam(klon, klev), knyam(klon, klev) !----------------------------------------------------------------------- - if (.not. (iflag_pbl >= 6 .and. iflag_pbl <= 9)) then - print *, 'probleme de coherence dans appel a MY' - stop 1 - endif + call assert(iflag_pbl >= 6 .and. iflag_pbl <= 9, "yamada4") ipas = ipas+1 ! les increments verticaux DO ig = 1, ngrid ! alerte: zlev n'est pas declare a nlev - 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)) ENDDO - DO k = 1, nlay + DO k = 1, klev DO ig = 1, ngrid unsdz(ig, k) = 1.E+0/(zlev(ig, k+1)-zlev(ig, k)) ENDDO ENDDO + DO ig = 1, ngrid unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1)) ENDDO - DO k = 2, nlay + + DO k = 2, klev DO ig = 1, ngrid unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1)) ENDDO ENDDO + DO ig = 1, ngrid - 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)) ENDDO do k = 2, klev @@ -132,7 +130,7 @@ else rif(ig, k) = rifc endif - if(rif(ig, k).lt.0.16) then + if (rif(ig, k).lt.0.16) then alpha(ig, k) = falpha(rif(ig, k)) sm(ig, k) = fsm(rif(ig, k)) else @@ -143,8 +141,8 @@ enddo enddo - ! Au premier appel, on détermine l et q2 de façon itérative. - ! Itération pour déterminer la longueur de mélange + ! Au premier appel, on d\'etermine l et q2 de façon it\'erative. + ! It\'eration pour d\'eterminer la longueur de m\'elange if (first .or. iflag_pbl == 6) then do ig = 1, ngrid @@ -199,7 +197,7 @@ do k = 2, klev do ig = 1, ngrid l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k)) - if(first) then + if (first) then q2(ig, k) = l(ig, k)**2 * zz(ig, k) endif enddo @@ -265,7 +263,6 @@ delta(ig, k) = 1.e-20 endif km(ig, k) = l(ig, k)*sqrt(q2(ig, k))*sm(ig, k) - aa0 = (m2(ig, k)-alpha(ig, k)*n2(ig, k)-delta(ig, k)/b1) aa1 = (m2(ig, k)*(1.-rif(ig, k))-delta(ig, k)/b1) aa(ig, k) = aa1*dt/(delta(ig, k)*l(ig, k)) qpre = sqrt(q2(ig, k)) @@ -288,7 +285,7 @@ enddo endif - ! Calcul des coefficients de mélange + ! Calcul des coefficients de m\'elange do k = 2, klev do ig = 1, ngrid zq = sqrt(q2(ig, k)) @@ -311,7 +308,7 @@ print *, 'pblhmin ', pblhmin do k = 2, klev - do ig = 1, klon + do ig = 1, ngrid 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 @@ -329,23 +326,14 @@ enddo enddo - ! 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) - first = .false. end SUBROUTINE yamada4 !******************************************************************* - function frif(ri) + real function frif(ri) - real frif real, intent(in):: ri frif = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156)) @@ -354,9 +342,8 @@ !******************************************************************* - function falpha(ri) + real function falpha(ri) - real falpha real, intent(in):: ri falpha = 1.318*(0.2231-ri)/(0.2341-ri) @@ -365,9 +352,8 @@ !******************************************************************* - function fsm(ri) + real function fsm(ri) - real fsm real, intent(in):: ri fsm = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri)) @@ -376,9 +362,8 @@ !******************************************************************* - function fl(zzz, zl0, zq2, zn2) + real function fl(zzz, zl0, zq2, zn2) - real fl real, intent(in):: zzz, zl0, zq2, zn2 fl = max(min(zl0 * kap * zzz / (kap * zzz + zl0), &