! $Header$ SUBROUTINE advn(q, masse, w, pbaru, pbarv, pdt, mode) ! ! Auteur : F. Hourdin ! ! ******************************************************************** ! Shema d'advection " pseudo amont " . ! ******************************************************************** ! q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... ! ! pbaru,pbarv,w flux de masse en u ,v ,w ! pdt pas de temps ! ! -------------------------------------------------------------------- USE lmdz_iniprint, ONLY: lunout, prt_level IMPLICIT NONE ! include "dimensions.h" include "paramet.h" include "comgeom.h" ! ! Arguments: ! ---------- INTEGER :: mode REAL :: masse(ip1jmp1, llm) REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) REAL :: q(ip1jmp1, llm) REAL :: w(ip1jmp1, llm), pdt ! ! Local ! --------- ! INTEGER :: i, ij, l, j, ii INTEGER :: ijlqmin, iqmin, jqmin, lqmin INTEGER :: ismin ! REAL :: zm(ip1jmp1, llm), newmasse REAL :: mu(ip1jmp1, llm) REAL :: mv(ip1jm, llm) REAL :: mw(ip1jmp1, llm + 1) REAL :: zq(ip1jmp1, llm), zz, qpn, qps REAL :: zqg(ip1jmp1, llm), zqd(ip1jmp1, llm) REAL :: zqs(ip1jmp1, llm), zqn(ip1jmp1, llm) REAL :: zqh(ip1jmp1, llm), zqb(ip1jmp1, llm) REAL :: temps0, temps1, temps2, temps3 REAL :: ztemps1, ztemps2, ssum save temps1, temps2, temps3 REAL :: zzpbar, zzw REAL :: qmin, qmax data qmin, qmax/0., 1./ data temps1, temps2, temps3/0., 0., 0./ zzpbar = 0.5 * pdt zzw = pdt DO l = 1, llm DO ij = iip2, ip1jm mu(ij, l) = pbaru(ij, l) * zzpbar ENDDO DO ij = 1, ip1jm mv(ij, l) = pbarv(ij, l) * zzpbar ENDDO DO ij = 1, ip1jmp1 mw(ij, l) = w(ij, l) * zzw ENDDO ENDDO DO ij = 1, ip1jmp1 mw(ij, llm + 1) = 0. ENDDO do l = 1, llm qpn = 0. qps = 0. do ij = 1, iim qpn = qpn + q(ij, l) * masse(ij, l) qps = qps + q(ip1jm + ij, l) * masse(ip1jm + ij, l) enddo qpn = qpn / ssum(iim, masse(1, l), 1) qps = qps / ssum(iim, masse(ip1jm + 1, l), 1) do ij = 1, iip1 q(ij, l) = qpn q(ip1jm + ij, l) = qps enddo enddo do ij = 1, ip1jmp1 mw(ij, llm + 1) = 0. enddo do l = 1, llm do ij = 1, ip1jmp1 zq(ij, l) = q(ij, l) zm(ij, l) = masse(ij, l) enddo enddo ! CALL minmaxq(zq,qmin,qmax,'avant vlx ') CALL advnqx(zq, zqg, zqd) CALL advnx(zq, zqg, zqd, zm, mu, mode) CALL advnqy(zq, zqs, zqn) CALL advny(zq, zqs, zqn, zm, mv) CALL advnqz(zq, zqh, zqb) CALL advnz(zq, zqh, zqb, zm, mw) ! CALL vlz(zq,0.,zm,mw) CALL advnqy(zq, zqs, zqn) CALL advny(zq, zqs, zqn, zm, mv) CALL advnqx(zq, zqg, zqd) CALL advnx(zq, zqg, zqd, zm, mu, mode) ! CALL minmaxq(zq,qmin,qmax,'apres vlx ') do l = 1, llm do ij = 1, ip1jmp1 q(ij, l) = zq(ij, l) enddo do ij = 1, ip1jm + 1, iip1 q(ij + iim, l) = q(ij, l) enddo enddo RETURN END SUBROUTINE advn SUBROUTINE advnqx(q, qg, qd) ! ! Auteurs: Calcul des valeurs de q aux point u. ! ! -------------------------------------------------------------------- USE lmdz_iniprint, ONLY: lunout, prt_level IMPLICIT NONE ! INCLUDE "dimensions.h" INCLUDE "paramet.h" ! ! ! Arguments: ! ---------- REAL :: q(ip1jmp1, llm), qg(ip1jmp1, llm), qd(ip1jmp1, llm) ! ! Local ! --------- ! INTEGER :: ij, l ! REAL :: dxqu(ip1jmp1), zqu(ip1jmp1) REAL :: zqmax(ip1jmp1), zqmin(ip1jmp1) LOGICAL :: extremum(ip1jmp1) INTEGER :: mode save mode data mode/1/ ! calcul des pentes en u: ! ----------------------- IF (mode==0) THEN do l = 1, llm do ij = 1, ip1jm qd(ij, l) = q(ij, l) qg(ij, l) = q(ij, l) enddo enddo else do l = 1, llm do ij = iip2, ip1jm - 1 dxqu(ij) = q(ij + 1, l) - q(ij, l) zqu(ij) = 0.5 * (q(ij + 1, l) + q(ij, l)) enddo do ij = iip1 + iip1, ip1jm, iip1 dxqu(ij) = dxqu(ij - iim) zqu(ij) = zqu(ij - iim) enddo do ij = iip2, ip1jm - 1 zqu(ij) = zqu(ij) - dxqu(ij + 1) / 12. enddo do ij = iip1 + iip1, ip1jm, iip1 zqu(ij) = zqu(ij - iim) enddo do ij = iip2 + 1, ip1jm zqu(ij) = zqu(ij) + dxqu(ij - 1) / 12. enddo do ij = iip1 + iip1, ip1jm, iip1 zqu(ij - iim) = zqu(ij) enddo ! calcul des valeurs max et min acceptees aux interfaces do ij = iip2, ip1jm - 1 zqmax(ij) = max(q(ij + 1, l), q(ij, l)) zqmin(ij) = min(q(ij + 1, l), q(ij, l)) enddo do ij = iip1 + iip1, ip1jm, iip1 zqmax(ij) = zqmax(ij - iim) zqmin(ij) = zqmin(ij - iim) enddo do ij = iip2 + 1, ip1jm extremum(ij) = dxqu(ij) * dxqu(ij - 1)<=0. enddo do ij = iip1 + iip1, ip1jm, iip1 extremum(ij - iim) = extremum(ij) enddo do ij = iip2, ip1jm zqu(ij) = min(max(zqmin(ij), zqu(ij)), zqmax(ij)) enddo do ij = iip2 + 1, ip1jm IF(extremum(ij)) THEN qg(ij, l) = q(ij, l) qd(ij, l) = q(ij, l) else qd(ij, l) = zqu(ij) qg(ij, l) = zqu(ij - 1) endif enddo do ij = iip1 + iip1, ip1jm, iip1 qd(ij - iim, l) = qd(ij, l) qg(ij - iim, l) = qg(ij, l) enddo goto 8888 do ij = iip2 + 1, ip1jm IF(extremum(ij).and..not.extremum(ij - 1)) & qd(ij - 1, l) = q(ij, l) enddo do ij = iip1 + iip1, ip1jm, iip1 qd(ij - iim, l) = qd(ij, l) enddo do ij = iip2, ip1jm - 1 IF (extremum(ij).and..not.extremum(ij + 1)) & qg(ij + 1, l) = q(ij, l) enddo do ij = iip1 + iip1, ip1jm, iip1 qg(ij, l) = qg(ij - iim, l) enddo 8888 continue enddo ENDIF RETURN END SUBROUTINE advnqx SUBROUTINE advnqy(q, qs, qn) ! ! Auteurs: Calcul des valeurs de q aux point v. ! ! -------------------------------------------------------------------- USE lmdz_iniprint, ONLY: lunout, prt_level IMPLICIT NONE ! INCLUDE "dimensions.h" INCLUDE "paramet.h" ! ! ! Arguments: ! ---------- REAL :: q(ip1jmp1, llm), qs(ip1jmp1, llm), qn(ip1jmp1, llm) ! ! Local ! --------- ! INTEGER :: ij, l ! REAL :: dyqv(ip1jm), zqv(ip1jm, llm) REAL :: zqmax(ip1jm), zqmin(ip1jm) LOGICAL :: extremum(ip1jmp1) INTEGER :: mode save mode data mode/1/ IF (mode==0) THEN do l = 1, llm do ij = 1, ip1jmp1 qn(ij, l) = q(ij, l) qs(ij, l) = q(ij, l) enddo enddo else ! calcul des pentes en u: ! ----------------------- do l = 1, llm do ij = 1, ip1jm dyqv(ij) = q(ij, l) - q(ij + iip1, l) enddo do ij = iip2, ip1jm - iip1 zqv(ij, l) = 0.5 * (q(ij + iip1, l) + q(ij, l)) zqv(ij, l) = zqv(ij, l) + (dyqv(ij + iip1) - dyqv(ij - iip1)) / 12. enddo do ij = iip2, ip1jm extremum(ij) = dyqv(ij) * dyqv(ij - iip1)<=0. enddo ! Pas de pentes aux poles do ij = 1, iip1 zqv(ij, l) = q(ij, l) zqv(ip1jm - iip1 + ij, l) = q(ip1jm + ij, l) extremum(ij) = .TRUE. extremum(ip1jmp1 - iip1 + ij) = .TRUE. enddo ! calcul des valeurs max et min acceptees aux interfaces do ij = 1, ip1jm zqmax(ij) = max(q(ij + iip1, l), q(ij, l)) zqmin(ij) = min(q(ij + iip1, l), q(ij, l)) enddo do ij = 1, ip1jm zqv(ij, l) = min(max(zqmin(ij), zqv(ij, l)), zqmax(ij)) enddo do ij = iip2, ip1jm IF(extremum(ij)) THEN qs(ij, l) = q(ij, l) qn(ij, l) = q(ij, l) ! if (.NOT.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l) ! if (.NOT.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l) else qs(ij, l) = zqv(ij, l) qn(ij, l) = zqv(ij - iip1, l) endif enddo do ij = 1, iip1 qs(ij, l) = q(ij, l) qn(ij, l) = q(ij, l) qs(ip1jm + ij, l) = q(ip1jm + ij, l) qn(ip1jm + ij, l) = q(ip1jm + ij, l) enddo enddo ENDIF RETURN END SUBROUTINE advnqy SUBROUTINE advnqz(q, qh, qb) ! ! Auteurs: Calcul des valeurs de q aux point v. ! ! -------------------------------------------------------------------- USE lmdz_iniprint, ONLY: lunout, prt_level IMPLICIT NONE ! INCLUDE "dimensions.h" INCLUDE "paramet.h" ! ! ! Arguments: ! ---------- REAL :: q(ip1jmp1, llm), qh(ip1jmp1, llm), qb(ip1jmp1, llm) ! ! Local ! --------- ! INTEGER :: ij, l ! REAL :: dzqw(ip1jmp1, llm + 1), zqw(ip1jmp1, llm + 1) REAL :: zqmax(ip1jmp1, llm), zqmin(ip1jmp1, llm) LOGICAL :: extremum(ip1jmp1, llm) INTEGER :: mode save mode data mode/1/ ! calcul des pentes en u: ! ----------------------- IF (mode==0) THEN do l = 1, llm do ij = 1, ip1jmp1 qb(ij, l) = q(ij, l) qh(ij, l) = q(ij, l) enddo enddo else do l = 2, llm do ij = 1, ip1jmp1 dzqw(ij, l) = q(ij, l - 1) - q(ij, l) zqw(ij, l) = 0.5 * (q(ij, l - 1) + q(ij, l)) enddo enddo do ij = 1, ip1jmp1 dzqw(ij, 1) = 0. dzqw(ij, llm + 1) = 0. enddo do l = 2, llm do ij = 1, ip1jmp1 zqw(ij, l) = zqw(ij, l) + (dzqw(ij, l + 1) - dzqw(ij, l - 1)) / 12. enddo enddo do l = 2, llm - 1 do ij = 1, ip1jmp1 extremum(ij, l) = dzqw(ij, l) * dzqw(ij, l + 1)<=0. enddo enddo ! Pas de pentes en bas et en haut do ij = 1, ip1jmp1 zqw(ij, 2) = q(ij, 1) zqw(ij, llm) = q(ij, llm) extremum(ij, 1) = .TRUE. extremum(ij, llm) = .TRUE. enddo ! calcul des valeurs max et min acceptees aux interfaces do l = 2, llm do ij = 1, ip1jmp1 zqmax(ij, l) = max(q(ij, l - 1), q(ij, l)) zqmin(ij, l) = min(q(ij, l - 1), q(ij, l)) enddo enddo do l = 2, llm do ij = 1, ip1jmp1 zqw(ij, l) = min(max(zqmin(ij, l), zqw(ij, l)), zqmax(ij, l)) enddo enddo do l = 2, llm - 1 do ij = 1, ip1jmp1 IF(extremum(ij, l)) THEN qh(ij, l) = q(ij, l) qb(ij, l) = q(ij, l) else qh(ij, l) = zqw(ij, l + 1) qb(ij, l) = zqw(ij, l) endif enddo enddo ! do l=2,llm-1 ! do ij=1,ip1jmp1 ! IF(extremum(ij,l)) THEN ! if (.NOT.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l) ! if (.NOT.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l) ! endif ! enddo ! enddo do ij = 1, ip1jmp1 qb(ij, 1) = q(ij, 1) qh(ij, 1) = q(ij, 1) qb(ij, llm) = q(ij, llm) qh(ij, llm) = q(ij, llm) enddo ENDIF RETURN END SUBROUTINE advnqz SUBROUTINE advnx(q, qg, qd, masse, u_m, mode) ! ! Auteur : F. Hourdin ! ! ******************************************************************** ! Shema d'advection " pseudo amont " . ! ******************************************************************** ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... ! ! ! -------------------------------------------------------------------- USE lmdz_iniprint, ONLY: lunout, prt_level IMPLICIT NONE ! include "dimensions.h" include "paramet.h" ! ! ! Arguments: ! ---------- INTEGER :: mode REAL :: masse(ip1jmp1, llm) REAL :: u_m(ip1jmp1, llm) REAL :: q(ip1jmp1, llm), qd(ip1jmp1, llm), qg(ip1jmp1, llm) ! ! Local ! --------- ! INTEGER :: i, j, ij, l, indu(ip1jmp1), niju, iju, ijq INTEGER :: n0, nl(llm) ! REAL :: new_m, zu_m, zdq, zz REAL :: zsigg(ip1jmp1, llm), zsigd(ip1jmp1, llm), zsig REAL :: u_mq(ip1jmp1, llm) REAL :: zm, zq, zsigm, zsigp, zqm, zqp, zu LOGICAL :: ladvplus(ip1jmp1, llm) REAL :: prec save prec data prec/1.e-15/ do l = 1, llm do ij = iip2, ip1jm zdq = qd(ij, l) - qg(ij, l) ! if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) THEN ! PRINT*,'probleme au point ij=',ij,' l=',l ! PRINT*,qd(ij,l),q(ij,l),qg(ij,l) ! qd(ij,l)=q(ij,l) ! qg(ij,l)=q(ij,l) ! END IF IF(abs(zdq)>prec) THEN zsigd(ij, l) = (q(ij, l) - qg(ij, l)) / zdq zsigg(ij, l) = 1. - zsigd(ij, l) ! IF(.NOT.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .AND. ! s zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) THEN ! PRINT*,'probleme au point ij=',ij,' l=',l ! PRINT*,'sigg=',zsigg(ij,l),' sigd=',zsigd(ij,l) ! PRINT*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq ! stop ! END IF else zsigd(ij, l) = 0.5 zsigg(ij, l) = 0.5 qd(ij, l) = q(ij, l) qg(ij, l) = q(ij, l) endif enddo enddo ! calcul de la pente maximum dans la maille en valeur absolue do l = 1, llm do ij = iip2, ip1jm - 1 IF (u_m(ij, l)>=0.) THEN zsigp = zsigd(ij, l) zsigm = zsigg(ij, l) zqp = qd(ij, l) zqm = qg(ij, l) zm = masse(ij, l) zq = q(ij, l) else zsigm = zsigd(ij + 1, l) zsigp = zsigg(ij + 1, l) zqm = qd(ij + 1, l) zqp = qg(ij + 1, l) zm = masse(ij + 1, l) zq = q(ij + 1, l) endif zu = abs(u_m(ij, l)) ladvplus(ij, l) = zu>zm zsig = zu / zm IF(zsig==0.) zsigp = 0.1 IF (mode==1) THEN IF (zsig<=zsigp) THEN u_mq(ij, l) = u_m(ij, l) * zqp ELSE IF (mode==1) THEN u_mq(ij, l) = & sign(zm, u_m(ij, l)) * (zsigp * zqp + (zsig - zsigp) * zqm) endif else IF (zsig<=zsigp) THEN u_mq(ij, l) = u_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq)) else zz = 0.5 * (zsig - zsigp) / zsigm u_mq(ij, l) = sign(zm, u_m(ij, l)) * (0.5 * (zq + zqp) * zsigp & + (zsig - zsigp) * (zq + zz * (zqm - zq))) endif endif ! IF(zsig.lt.0.) THEN ! PRINT*,'au point ij=',ij,' l=',l,' sig=',zsig ! stop ! END IF enddo enddo do l = 1, llm do ij = iip1 + iip1, ip1jm, iip1 u_mq(ij, l) = u_mq(ij - iim, l) ladvplus(ij, l) = ladvplus(ij - iim, l) enddo enddo !================================================================= ! SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES !================================================================= ! tris des regions a traiter n0 = 0 do l = 1, llm nl(l) = 0 do ij = iip2, ip1jm IF(ladvplus(ij, l)) THEN nl(l) = nl(l) + 1 u_mq(ij, l) = 0. endif enddo n0 = n0 + nl(l) enddo IF(n0>1) THEN IF (prt_level > 9) WRITE(lunout, *) & 'Nombre de points pour lesquels on advect plus que le' & , 'contenu de la maille : ', n0 do l = 1, llm IF(nl(l)>0) THEN iju = 0 ! indicage des mailles concernees par le traitement special do ij = iip2, ip1jm IF(ladvplus(ij, l).AND.mod(ij, iip1)/=0) THEN iju = iju + 1 indu(iju) = ij endif enddo niju = iju ! PRINT*,'niju,nl',niju,nl(l) ! traitement des mailles do iju = 1, niju ij = indu(iju) j = (ij - 1) / iip1 + 1 zu_m = u_m(ij, l) u_mq(ij, l) = 0. IF(zu_m>0.) THEN ijq = ij i = ijq - (j - 1) * iip1 ! accumulation pour les mailles completements advectees do while(zu_m>masse(ijq, l)) u_mq(ij, l) = u_mq(ij, l) + q(ijq, l) * masse(ijq, l) zu_m = zu_m - masse(ijq, l) i = mod(i - 2 + iim, iim) + 1 ijq = (j - 1) * iip1 + i enddo ! MODIFS SPECIFIQUES DU SCHEMA ! ajout de la maille non completement advectee zsig = zu_m / masse(ijq, l) IF(zsig<=zsigd(ijq, l)) THEN u_mq(ij, l) = u_mq(ij, l) + zu_m * (qd(ijq, l) & - 0.5 * zsig / zsigd(ijq, l) * (qd(ijq, l) - q(ijq, l))) else ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l) ! goto 8888 zz = 0.5 * (zsig - zsigd(ijq, l)) / zsigg(ijq, l) IF(.NOT.(zz>0..and.zz<=0.5)) THEN WRITE(lunout, *)'probleme2 au point ij=', ij, & ' l=', l WRITE(lunout, *)'zz=', zz stop endif u_mq(ij, l) = u_mq(ij, l) + masse(ijq, l) * (& 0.5 * (q(ijq, l) + qd(ijq, l)) * zsigd(ijq, l) & + (zsig - zsigd(ijq, l)) * (q(ijq, l) + zz * (qg(ijq, l) - q(ijq, l)))) endif else ijq = ij + 1 i = ijq - (j - 1) * iip1 ! accumulation pour les mailles completements advectees do while(-zu_m>masse(ijq, l)) u_mq(ij, l) = u_mq(ij, l) - q(ijq, l) * masse(ijq, l) zu_m = zu_m + masse(ijq, l) i = mod(i, iim) + 1 ijq = (j - 1) * iip1 + i enddo ! ajout de la maille non completement advectee ! 2eme MODIF SPECIFIQUE zsig = -zu_m / masse(ij + 1, l) IF(zsig<=zsigg(ijq, l)) THEN u_mq(ij, l) = u_mq(ij, l) + zu_m * (qg(ijq, l) & - 0.5 * zsig / zsigg(ijq, l) * (qg(ijq, l) - q(ijq, l))) else ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l) ! goto 9999 zz = 0.5 * (zsig - zsigg(ijq, l)) / zsigd(ijq, l) IF(.NOT.(zz>0..and.zz<=0.5)) THEN WRITE(lunout, *)'probleme22 au point ij=', ij & , ' l=', l WRITE(lunout, *)'zz=', zz stop endif u_mq(ij, l) = u_mq(ij, l) - masse(ijq, l) * (& 0.5 * (q(ijq, l) + qg(ijq, l)) * zsigg(ijq, l) & + (zsig - zsigg(ijq, l)) * & (q(ijq, l) + zz * (qd(ijq, l) - q(ijq, l)))) endif ! fin de la modif endif enddo endif enddo ENDIF ! n0.gt.0 ! bouclage en latitude do l = 1, llm do ij = iip1 + iip1, ip1jm, iip1 u_mq(ij, l) = u_mq(ij - iim, l) enddo enddo !================================================================= ! CALCUL DE LA CONVERGENCE DES FLUX !================================================================= do l = 1, llm do ij = iip2 + 1, ip1jm new_m = masse(ij, l) + u_m(ij - 1, l) - u_m(ij, l) q(ij, l) = (q(ij, l) * masse(ij, l) + & u_mq(ij - 1, l) - u_mq(ij, l)) & / new_m masse(ij, l) = new_m enddo ! Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous) do ij = iip1 + iip1, ip1jm, iip1 q(ij - iim, l) = q(ij, l) masse(ij - iim, l) = masse(ij, l) enddo enddo RETURN END SUBROUTINE advnx SUBROUTINE advny(q, qs, qn, masse, v_m) ! ! Auteur : F. Hourdin ! ! ******************************************************************** ! Shema d'advection " pseudo amont " . ! ******************************************************************** ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... ! ! ! -------------------------------------------------------------------- USE lmdz_iniprint, ONLY: lunout, prt_level IMPLICIT NONE ! INCLUDE "dimensions.h" INCLUDE "paramet.h" INCLUDE "comgeom.h" ! ! ! Arguments: ! ---------- REAL :: masse(ip1jmp1, llm) REAL :: v_m(ip1jm, llm) REAL :: q(ip1jmp1, llm), qn(ip1jmp1, llm), qs(ip1jmp1, llm) ! ! Local ! --------- ! INTEGER :: ij, l ! REAL :: new_m, zdq, zz REAL :: zsigs(ip1jmp1), zsign(ip1jmp1), zsig REAL :: v_mq(ip1jm, llm) REAL :: convpn, convps, convmpn, convmps, massen, masses REAL :: zm, zq, zsigm, zsigp, zqm, zqp REAL :: ssum REAL :: prec save prec data prec/1.e-15/ do l = 1, llm do ij = 1, ip1jmp1 zdq = qn(ij, l) - qs(ij, l) ! if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) THEN ! PRINT*,'probleme au point ij=',ij,' l=',l,' advnqx' ! PRINT*,qn(ij,l),q(ij,l),qs(ij,l) ! qn(ij,l)=q(ij,l) ! qs(ij,l)=q(ij,l) ! END IF IF(abs(zdq)>prec) THEN zsign(ij) = (q(ij, l) - qs(ij, l)) / zdq zsigs(ij) = 1. - zsign(ij) ! IF(.NOT.(zsign(ij).ge.0..and.zsign(ij).le.1. .AND. ! s zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) THEN ! PRINT*,'probleme au point ij=',ij,' l=',l ! PRINT*,'sigs=',zsigs(ij),' sign=',zsign(ij) ! stop ! END IF else zsign(ij) = 0.5 zsigs(ij) = 0.5 endif enddo ! calcul de la pente maximum dans la maille en valeur absolue do ij = 1, ip1jm IF (v_m(ij, l)>=0.) THEN zsigp = zsign(ij + iip1) zsigm = zsigs(ij + iip1) zqp = qn(ij + iip1, l) zqm = qs(ij + iip1, l) zm = masse(ij + iip1, l) zq = q(ij + iip1, l) else zsigm = zsign(ij) zsigp = zsigs(ij) zqm = qn(ij, l) zqp = qs(ij, l) zm = masse(ij, l) zq = q(ij, l) endif zsig = abs(v_m(ij, l)) / zm IF(zsig==0.) zsigp = 0.1 IF (zsig<=zsigp) THEN v_mq(ij, l) = v_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq)) else zz = 0.5 * (zsig - zsigp) / zsigm v_mq(ij, l) = sign(zm, v_m(ij, l)) * (0.5 * (zq + zqp) * zsigp & + (zsig - zsigp) * (zq + zz * (zqm - zq))) endif enddo enddo do l = 1, llm do ij = iip2, ip1jm new_m = masse(ij, l) & + v_m(ij, l) - v_m(ij - iip1, l) q(ij, l) = (q(ij, l) * masse(ij, l) + v_mq(ij, l) - v_mq(ij - iip1, l)) & / new_m masse(ij, l) = new_m enddo !.-. ancienne version convpn = SSUM(iim, v_mq(1, l), 1) convmpn = ssum(iim, v_m(1, l), 1) massen = ssum(iim, masse(1, l), 1) new_m = massen + convmpn q(1, l) = (q(1, l) * massen + convpn) / new_m do ij = 1, iip1 q(ij, l) = q(1, l) masse(ij, l) = new_m * aire(ij) / apoln enddo convps = -SSUM(iim, v_mq(ip1jm - iim, l), 1) convmps = -ssum(iim, v_m(ip1jm - iim, l), 1) masses = ssum(iim, masse(ip1jm + 1, l), 1) new_m = masses + convmps q(ip1jm + 1, l) = (q(ip1jm + 1, l) * masses + convps) / new_m do ij = ip1jm + 1, ip1jmp1 q(ij, l) = q(ip1jm + 1, l) masse(ij, l) = new_m * aire(ij) / apols enddo enddo RETURN END SUBROUTINE advny SUBROUTINE advnz(q, qh, qb, masse, w_m) ! ! Auteurs: F.Hourdin ! ! ******************************************************************** ! Shema d'advection " pseudo amont " . ! b designe le bas et h le haut ! il y a une correspondance entre le b en z et le d en x ! ******************************************************************** ! ! ! -------------------------------------------------------------------- USE lmdz_iniprint, ONLY: lunout, prt_level IMPLICIT NONE ! INCLUDE "dimensions.h" INCLUDE "paramet.h" INCLUDE "comgeom.h" ! ! ! Arguments: ! ---------- REAL :: masse(ip1jmp1, llm) REAL :: w_m(ip1jmp1, llm + 1) REAL :: q(ip1jmp1, llm), qb(ip1jmp1, llm), qh(ip1jmp1, llm) ! ! Local ! --------- ! INTEGER :: ij, l ! REAL :: new_m, zdq, zz REAL :: zsigh(ip1jmp1, llm), zsigb(ip1jmp1, llm), zsig REAL :: w_mq(ip1jmp1, llm + 1) REAL :: zm, zq, zsigm, zsigp, zqm, zqp REAL :: prec save prec data prec/1.e-13/ do l = 1, llm do ij = 1, ip1jmp1 zdq = qb(ij, l) - qh(ij, l) ! if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) THEN ! PRINT*,'probleme au point ij=',ij,' l=',l ! PRINT*,qh(ij,l),q(ij,l),qb(ij,l) ! qh(ij,l)=q(ij,l) ! qb(ij,l)=q(ij,l) ! END IF IF(abs(zdq)>prec) THEN zsigb(ij, l) = (q(ij, l) - qh(ij, l)) / zdq zsigh(ij, l) = 1. - zsigb(ij, l) zsigb(ij, l) = min(max(zsigb(ij, l), 0.), 1.) else zsigb(ij, l) = 0.5 zsigh(ij, l) = 0.5 endif enddo enddo ! PRINT*,'ok1' ! calcul de la pente maximum dans la maille en valeur absolue do l = 2, llm do ij = 1, ip1jmp1 IF (w_m(ij, l)>=0.) THEN zsigp = zsigb(ij, l) zsigm = zsigh(ij, l) zqp = qb(ij, l) zqm = qh(ij, l) zm = masse(ij, l) zq = q(ij, l) else zsigm = zsigb(ij, l - 1) zsigp = zsigh(ij, l - 1) zqm = qb(ij, l - 1) zqp = qh(ij, l - 1) zm = masse(ij, l - 1) zq = q(ij, l - 1) endif zsig = abs(w_m(ij, l)) / zm IF(zsig==0.) zsigp = 0.1 IF (zsig<=zsigp) THEN w_mq(ij, l) = w_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq)) else zz = 0.5 * (zsig - zsigp) / zsigm w_mq(ij, l) = sign(zm, w_m(ij, l)) * (0.5 * (zq + zqp) * zsigp & + (zsig - zsigp) * (zq + zz * (zqm - zq))) endif enddo enddo do ij = 1, ip1jmp1 w_mq(ij, llm + 1) = 0. w_mq(ij, 1) = 0. enddo do l = 1, llm do ij = 1, ip1jmp1 new_m = masse(ij, l) + w_m(ij, l + 1) - w_m(ij, l) q(ij, l) = (q(ij, l) * masse(ij, l) + w_mq(ij, l + 1) - w_mq(ij, l)) & / new_m masse(ij, l) = new_m enddo enddo ! PRINT*,'ok3' RETURN END SUBROUTINE advnz