--- trunk/libf/phylmd/coefkz.f90 2012/04/20 14:58:43 61 +++ trunk/libf/phylmd/coefkz.f90 2012/07/26 14:37:37 62 @@ -5,20 +5,21 @@ contains SUBROUTINE coefkz(nsrf, knon, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, & - t, q, qsurf, pcfm, pcfh) + t, q, qsurf, coefm, coefh) - ! Authors: F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) + ! Authors: F. Hourdin, M. Forichon, Z. X. Li (LMD/CNRS) ! date: 1993/09/22 ! Objet : calculer le coefficient de frottement du sol ("Cdrag") et les ! coefficients d'échange turbulent dans l'atmosphère. - USE indicesol, ONLY : is_oce - USE dimphy, ONLY : klev, klon, max - USE conf_gcm_m, ONLY : prt_level - USE suphec_m, ONLY : rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt - USE yoethf_m, ONLY : r2es, r5ies, r5les, rvtmp2 - USE fcttre, ONLY : dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep - USE conf_phys_m, ONLY : iflag_pbl + USE indicesol, ONLY: is_oce + USE dimphy, ONLY: klev, klon, max + USE conf_gcm_m, ONLY: prt_level + USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt + USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2 + USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep + USE conf_phys_m, ONLY: iflag_pbl + use clcdrag_m, only: clcdrag ! Arguments: @@ -38,12 +39,14 @@ REAL, intent(in):: t(klon, klev) ! temperature (K) real, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg) real, intent(in):: qsurf(klon) - REAL, intent(inout):: pcfm(klon, klev) ! coefficient, vitesse - real, intent(inout):: pcfh(klon, klev) ! coefficient, chaleur et humidité + REAL, intent(out):: coefm(:, :) ! (knon, klev) coefficient, vitesse + + real, intent(out):: coefh(:, :) ! (knon, klev) + ! coefficient, chaleur et humidité ! Local: - INTEGER itop(klon) ! numero de couche du sommet de la couche limite + INTEGER itop(knon) ! numero de couche du sommet de la couche limite ! Quelques constantes et options: @@ -53,14 +56,17 @@ REAL, PARAMETER:: cc = 5. REAL, PARAMETER:: cd = 5. REAL, PARAMETER:: clam = 160. - REAL, PARAMETER:: ratqs = 0.05 ! ! largeur de distribution de vapeur d'eau - LOGICAL, PARAMETER:: richum = .TRUE. ! utilise le nombre de Richardson humide + REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau + + LOGICAL, PARAMETER:: richum = .TRUE. + ! utilise le nombre de Richardson humide + REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique REAL, PARAMETER:: prandtl = 0.4 REAL kstable ! diffusion minimale (situation stable) - REAL, PARAMETER:: mixlen = 35. ! constante controlant longueur de melange - INTEGER, PARAMETER:: isommet = klev ! le sommet de la couche limite + REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange + INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite LOGICAL, PARAMETER:: tvirtu = .TRUE. ! calculer Ri d'une maniere plus performante @@ -71,32 +77,20 @@ INTEGER i, k, kk REAL zgeop(klon, klev) REAL zmgeom(klon) - REAL zri(klon) - REAL zl2(klon) + REAL ri(klon) + REAL l2(klon) REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon) - REAL pcfm1(klon), pcfh1(klon) - REAL zdphi, zdu2, ztvd, ztvu, zcdn - REAL zscf + REAL zdphi, zdu2, ztvd, ztvu, cdn + REAL scf REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs - REAL z2geomf, zalh2, zalm2, zscfh, zscfm + REAL z2geomf, zalh2, alm2, zscfh, scfm REAL, PARAMETER:: t_coup = 273.15 REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre !-------------------------------------------------------------------- - ! Initialiser les sorties - DO k = 1, klev - DO i = 1, knon - pcfm(i, k) = 0. - pcfh(i, k) = 0. - ENDDO - ENDDO - DO i = 1, knon - itop(i) = 0 - ENDDO - ! Prescrire la valeur de contre-gradient if (iflag_pbl.eq.1) then DO k = 3, klev @@ -139,18 +133,11 @@ ENDDO CALL clcdrag(klon, knon, nsrf, .false., u1, v1, t1, q1, z1, ts, qsurf, & - rugos, pcfm1, pcfh1) - - DO i = 1, knon - pcfm(i, 1) = pcfm1(i) - pcfh(i, 1) = pcfh1(i) - ENDDO + rugos, coefm(:, 1), coefh(:, 1)) ! Calculer les coefficients turbulents dans l'atmosphere - DO i = 1, knon - itop(i) = isommet - ENDDO + itop = isommet loop_vertical: DO k = 2, isommet loop_horiz: DO i = 1, knon @@ -173,7 +160,7 @@ zqs = zqs*zcor zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor) ELSE - IF (zt .LT. t_coup) THEN + IF (zt < t_coup) THEN zqs = qsats(zt) / pplay(i, k) zdqs = dqsats(zt, zqs) ELSE @@ -199,18 +186,18 @@ - zdphi/RCPD/(1.+RVTMP2*zq) & *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) & )*(1.+RETV*q(i, k-1)) - zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu)) - zri(i) = zri(i) & + ri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu)) + ri(i) = ri(i) & + zmgeom(i)*zmgeom(i)/RG*gamt(k) & *(paprs(i, k)/101325.0)**RKAPPA & /(zdu2*0.5*(ztvd+ztvu)) ELSE ! calcul de Ridchardson compatible LMD5 - zri(i) =(RCPD*(t(i, k)-t(i, k-1)) & + ri(i) =(RCPD*(t(i, k)-t(i, k-1)) & -RD*0.5*(t(i, k)+t(i, k-1))/paprs(i, k) & *(pplay(i, k)-pplay(i, k-1)) & )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i, k-1)+t(i, k))) - zri(i) = zri(i) + & + ri(i) = ri(i) + & zmgeom(i)*zmgeom(i)*gamt(k)/RG & *(paprs(i, k)/101325.0)**RKAPPA & /(zdu2*0.5*(t(i, k-1)+t(i, k))) @@ -218,49 +205,44 @@ ! finalement, les coefficients d'echange sont obtenus: - zcdn = SQRT(zdu2) / zmgeom(i) * RG + cdn = SQRT(zdu2) / zmgeom(i) * RG IF (opt_ec) THEN z2geomf = zgeop(i, k-1)+zgeop(i, k) - zalm2 = (0.5*ckap/RG*z2geomf & + alm2 = (0.5*ckap/RG*z2geomf & /(1.+0.5*ckap/rg/clam*z2geomf))**2 zalh2 = (0.5*ckap/rg*z2geomf & /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2 - IF (zri(i).LT.0.0) THEN + IF (ri(i) < 0.) THEN ! situation instable - zscf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 & + scf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 & / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG) - zscf = SQRT(-zri(i)*zscf) - zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf) - zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf) - pcfm(i, k) = zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm) - pcfh(i, k) = zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh) + scf = SQRT(-ri(i)*scf) + scfm = 1.0 / (1.0+3.0*cb*cc*alm2*scf) + zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*scf) + coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm) + coefh(i, k) = cdn*zalh2*(1.-3.0*cb*ri(i)*zscfh) ELSE ! situation stable - zscf = SQRT(1.+cd*zri(i)) - pcfm(i, k) = zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf) - pcfh(i, k) = zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf) + scf = SQRT(1.+cd*ri(i)) + coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf) + coefh(i, k) = cdn*zalh2/(1.+3.0*cb*ri(i)*scf) ENDIF ELSE - zl2(i) = (mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) & + l2(i) = (mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) & /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2 - pcfm(i, k) = sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable)) - pcfm(i, k)= zl2(i)* pcfm(i, k) - pcfh(i, k) = pcfm(i, k) /prandtl ! h et m different + coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable)) + coefm(i, k)= l2(i) * coefm(i, k) + coefh(i, k) = coefm(i, k) / prandtl ! h et m different ENDIF ENDDO loop_horiz ENDDO loop_vertical - ! Au-dela du sommet, pas de diffusion turbulente: - - DO i = 1, knon - IF (itop(i)+1 .LE. klev) THEN - DO k = itop(i)+1, klev - pcfh(i, k) = 0. - pcfm(i, k) = 0. - ENDDO - ENDIF - ENDDO + ! Au-delà du sommet, pas de diffusion turbulente : + forall (i = 1: knon) + coefh(i, itop(i) + 1:) = 0. + coefm(i, itop(i) + 1:) = 0. + END forall END SUBROUTINE coefkz