--- trunk/libf/dyn3d/Read_reanalyse/reanalyse2nat.f90 2010/12/21 15:45:48 37 +++ trunk/Sources/dyn3d/Guide/Read_reanalyse/reanalyse2nat.f 2015/10/06 15:57:02 173 @@ -1,146 +1,117 @@ +module reanalyse2nat_m + implicit none -!=========================================================================== - subroutine reanalyse2nat(nlevnc,psi & - ,unc,vnc,tnc,qnc,psnc,pl,u,v,t,q & - ,ps,masse,pk) -!=========================================================================== - -! ----------------------------------------------------------------- -! Inversion Nord/sud de la grille + interpollation sur les niveaux -! verticaux du modele. -! ----------------------------------------------------------------- - - use dimens_m - use paramet_m - use comconst - use comvert - use comgeom - use exner_hyb_m, only: exner_hyb - use guide_m - - implicit none - - - integer nlevnc - real psi(iip1,jjp1) - real u(iip1,jjp1,llm),v(iip1,jjm,llm) - real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm) - - real pl(nlevnc) - real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc) - real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1) - real qnc(iip1,jjp1,nlevnc) - - real zu(iip1,jjp1,llm),zv(iip1,jjm,llm) - real zt(iip1,jjp1,llm),zq(iip1,jjp1,llm) - - real pext(iip1,jjp1,llm) - real pbarx(iip1,jjp1,llm),pbary(iip1,jjm,llm) - real plunc(iip1,jjp1,llm),plvnc(iip1,jjm,llm) - real plsnc(iip1,jjp1,llm) - - real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm),pks(iip1,jjp1) - real pkf(iip1,jjp1,llm) - real masse(iip1,jjp1,llm),pls(iip1,jjp1,llm) - real prefkap,unskap - - - integer i,j,l - - -! ----------------------------------------------------------------- -! calcul de la pression au milieu des couches. -! ----------------------------------------------------------------- - - forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * psi - call massdair(p,masse) - CALL exner_hyb(psi,p,pks,pk,pkf) - -! .... Calcul de pls , pression au milieu des couches ,en Pascals - unskap=1./kappa - prefkap = preff ** kappa -! PRINT *,' Pref kappa unskap ',preff,kappa,unskap - DO l = 1, llm - DO j=1,jjp1 - DO i =1, iip1 - pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap - ENDDO - ENDDO +contains + + subroutine reanalyse2nat(invert_y, psi, unc, vnc, tnc, qnc, pl, u, v, t, q, & + pk) + + ! Inversion nord-sud de la grille et interpolation verticale sur + ! les niveaux du modèle. + + USE comconst, ONLY: cpp, kappa + USE comgeom, ONLY: aireu_2d, airev_2d, aire_2d + USE dimens_m, ONLY: iim, jjm, llm + USE disvert_m, ONLY: ap, bp, preff + USE exner_hyb_m, ONLY: exner_hyb + use massbar_m, only: massbar + use massdair_m, only: massdair + USE paramet_m, ONLY: iip1, jjp1, llmp1 + use pres2lev_m, only: pres2lev + + logical, intent(in):: invert_y + real, intent(in):: psi(:, :) ! (iip1, jjp1) + + real, intent(in):: unc(:, :, :) ! (iip1, jjp1, :) + real, intent(in):: vnc(:, :, :) ! (iip1, jjm, :) + real, intent(in):: tnc(:, :, :) ! (iip1, jjp1, :) + real, intent(in):: qnc(:, :, :) ! (iip1, jjp1, :) + real, intent(in):: pl(:) + + real, intent(out):: u(:, :, :) ! (iip1, jjp1, llm) + real, intent(out):: v(:, :, :) ! (iip1, jjm, llm) + real, intent(out):: t(:, :, :), q(:, :, :) ! (iip1, jjp1, llm) + real, intent(out):: pk(:, :, :) ! (iip1, jjp1, llm) + + ! Local: + + real zu(iip1, jjp1, llm), zv(iip1, jjm, llm) + real zt(iip1, jjp1, llm), zq(iip1, jjp1, llm) + + real pext(iip1, jjp1, llm) + real pbarx(iip1, jjp1, llm), pbary(iip1, jjm, llm) + real plunc(iip1, jjp1, llm), plvnc(iip1, jjm, llm) + real plsnc(iip1, jjp1, llm) + + real p(iip1, jjp1, llmp1) + real pks(iip1, jjp1) + real pls(iip1, jjp1, llm) + real prefkap, unskap + + integer i, j, l + + ! ----------------------------------------------------------------- + + ! calcul de la pression au milieu des couches + forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * psi + CALL exner_hyb(psi, p, pks, pk) + + ! Calcul de pls, pression au milieu des couches, en Pascals + unskap=1./kappa + prefkap = preff ** kappa + DO l = 1, llm + DO j=1, jjp1 + DO i =1, iip1 + pls(i, j, l) = preff * ( pk(i, j, l)/cpp) ** unskap + ENDDO ENDDO + ENDDO + + ! calcul des pressions pour les grilles u et v + do l=1, llm + do j=1, jjp1 + do i=1, iip1 + pext(i, j, l)=pls(i, j, l)*aire_2d(i, j) + enddo + enddo + enddo + call massbar(pext, pbarx, pbary ) + do l=1, llm + do j=1, jjp1 + do i=1, iip1 + plunc(i, jjp1+1-j, l)=pbarx(i, j, l)/aireu_2d(i, j) + plsnc(i, jjp1+1-j, l)=pls(i, j, l) + enddo + enddo + enddo + do l=1, llm + do j=1, jjm + do i=1, iip1 + plvnc(i, jjm+1-j, l)=pbary(i, j, l)/airev_2d(i, j) + enddo + enddo + enddo + + call pres2lev(unc, zu, pl, plunc) + call pres2lev(vnc, zv, pl, plvnc ) + call pres2lev(tnc, zt, pl, plsnc) + call pres2lev(qnc, zq, pl, plsnc) + + if (invert_y) then + ! Inversion Nord/Sud + u=zu(:, jjp1:1:-1, :) + v=zv(:, jjm:1:-1, :) + t=zt(:, jjp1:1:-1, :) + q=zq(:, jjp1:1:-1, :) + else + u = zu + v = zv + t = zt + q = zq + end if -! ----------------------------------------------------------------- -! calcul des pressions pour les grilles u et v -! ----------------------------------------------------------------- - - do l=1,llm - do j=1,jjp1 - do i=1,iip1 - pext(i,j,l)=pls(i,j,l)*aire_2d(i,j) - enddo - enddo - enddo - call massbar(pext, pbarx, pbary ) - do l=1,llm - do j=1,jjp1 - do i=1,iip1 - plunc(i,jjp1+1-j,l)=pbarx(i,j,l)/aireu_2d(i,j) - plsnc(i,jjp1+1-j,l)=pls(i,j,l) - enddo - enddo - enddo - do l=1,llm - do j=1,jjm - do i=1,iip1 - plvnc(i,jjm+1-j,l)=pbary(i,j,l)/airev_2d(i,j) - enddo - enddo - enddo - -! ----------------------------------------------------------------- - - if (guide_P) then - do j=1,jjp1 - do i=1,iim - ps(i,j)=psnc(i,jjp1+1-j) - enddo - ps(iip1,j)=ps(1,j) - enddo - endif - - -! ----------------------------------------------------------------- - call pres2lev(unc,zu,nlevnc,llm,pl,plunc,iip1,jjp1) - call pres2lev(vnc,zv,nlevnc,llm,pl,plvnc,iip1,jjm ) - call pres2lev(tnc,zt,nlevnc,llm,pl,plsnc,iip1,jjp1) - call pres2lev(qnc,zq,nlevnc,llm,pl,plsnc,iip1,jjp1) - -! call dump2d(iip1,jjp1,ps,'PS ') -! call dump2d(iip1,jjp1,psu,'PS ') -! call dump2d(iip1,jjm,psv,'PS ') -! Inversion Nord/Sud - do l=1,llm - do j=1,jjp1 - do i=1,iim - u(i,j,l)=zu(i,jjp1+1-j,l) - t(i,j,l)=zt(i,jjp1+1-j,l) - q(i,j,l)=zq(i,jjp1+1-j,l) - enddo - u(iip1,j,l)=u(1,j,l) - t(iip1,j,l)=t(1,j,l) - q(iip1,j,l)=q(1,j,l) - enddo - enddo - - do l=1,llm - do j=1,jjm - do i=1,iim - v(i,j,l)=zv(i,jjm+1-j,l) - enddo - v(iip1,j,l)=v(1,j,l) - enddo - enddo + end subroutine reanalyse2nat - return - end +end module reanalyse2nat_m