! $Id: vlsplt.F90 5159 2024-08-02 19:58:25Z abarral $ ! SUBROUTINE vlsplt(q, pente_max, masse, w, pbaru, pbarv, pdt, iq) USE infotrac, ONLY: nqtot, tracers USE lmdz_ssum_scopy, ONLY: scopy ! Auteurs: P.Le Van, F.Hourdin, F.Forget ! ******************************************************************** ! Shema d'advection " pseudo amont " . ! ******************************************************************** ! q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... ! pente_max facteur de limitation des pentes: 2 en general ! 0 pour un schema amont ! pbaru,pbarv,w flux de masse en u ,v ,w ! pdt pas de temps ! -------------------------------------------------------------------- USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE ! ! Arguments: ! ---------- REAL :: masse(ip1jmp1, llm), pente_max REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) REAL :: q(ip1jmp1, llm, nqtot) REAL :: w(ip1jmp1, llm), pdt INTEGER :: iq ! CRisi ! Local ! --------- INTEGER :: ij, l REAL :: zm(ip1jmp1, llm, nqtot) REAL :: mu(ip1jmp1, llm) REAL :: mv(ip1jm, llm) REAL :: mw(ip1jmp1, llm + 1) REAL :: zq(ip1jmp1, llm, nqtot) REAL :: zzpbar, zzw INTEGER :: ifils, iq2 ! CRisi REAL :: qmin, qmax DATA qmin, qmax/0., 1.e33/ 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 CALL SCOPY(ijp1llm, q(1, 1, iq), 1, zq(1, 1, iq), 1) CALL SCOPY(ijp1llm, masse, 1, zm(1, 1, iq), 1) DO ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) CALL SCOPY(ijp1llm, q(1, 1, iq2), 1, zq(1, 1, iq2), 1) enddo CALL vlx(zq, pente_max, zm, mu, iq) CALL vly(zq, pente_max, zm, mv, iq) CALL vlz(zq, pente_max, zm, mw, iq) CALL vly(zq, pente_max, zm, mv, iq) CALL vlx(zq, pente_max, zm, mu, iq) DO l = 1, llm DO ij = 1, ip1jmp1 q(ij, l, iq) = zq(ij, l, iq) ENDDO DO ij = 1, ip1jm + 1, iip1 q(ij + iim, l, iq) = q(ij, l, iq) ENDDO ENDDO ! CRisi: aussi pour les fils DO ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) DO l = 1, llm DO ij = 1, ip1jmp1 q(ij, l, iq2) = zq(ij, l, iq2) ENDDO DO ij = 1, ip1jm + 1, iip1 q(ij + iim, l, iq2) = q(ij, l, iq2) ENDDO ENDDO enddo END SUBROUTINE vlsplt RECURSIVE SUBROUTINE vlx(q, pente_max, masse, u_m, iq) USE infotrac, ONLY: nqtot, tracers, & ! CRisi min_qParent, min_qMass, min_ratio ! MVals et CRisi USE lmdz_iniprint, ONLY: lunout, prt_level ! Auteurs: P.Le Van, F.Hourdin, F.Forget ! ******************************************************************** ! Shema d'advection " pseudo amont " . ! ******************************************************************** ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... ! -------------------------------------------------------------------- USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE ! ! Arguments: ! ---------- REAL :: masse(ip1jmp1, llm, nqtot), pente_max REAL :: u_m(ip1jmp1, llm) REAL :: q(ip1jmp1, llm, nqtot) INTEGER :: iq ! CRisi ! Local ! --------- INTEGER :: ij, l, j, i, iju, ijq, indu(ip1jmp1), niju INTEGER :: n0, iadvplus(ip1jmp1, llm), nl(llm) REAL :: new_m, zu_m, zdum(ip1jmp1, llm) REAL :: dxq(ip1jmp1, llm), dxqu(ip1jmp1) REAL :: zz(ip1jmp1) REAL :: adxqu(ip1jmp1), dxqmax(ip1jmp1, llm) REAL :: u_mq(ip1jmp1, llm) ! CRisi REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) INTEGER :: ifils, iq2 ! CRisi LOGICAL, SAVE :: first DATA first/.TRUE./ ! calcul de la pente a droite et a gauche de la maille IF (pente_max>-1.e-5) THEN ! IF (pente_max.gt.10) THEN ! calcul des pentes avec limitation, Van Leer scheme I: ! ----------------------------------------------------- ! calcul de la pente aux points u DO l = 1, llm DO ij = iip2, ip1jm - 1 dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq) ENDDO DO ij = iip1 + iip1, ip1jm, iip1 dxqu(ij) = dxqu(ij - iim) ! sigu(ij)=sigu(ij-iim) ENDDO DO ij = iip2, ip1jm adxqu(ij) = abs(dxqu(ij)) ENDDO ! calcul de la pente maximum dans la maille en valeur absolue DO ij = iip2 + 1, ip1jm dxqmax(ij, l) = pente_max * & min(adxqu(ij - 1), adxqu(ij)) ! limitation subtile ! , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij))) ENDDO DO ij = iip1 + iip1, ip1jm, iip1 dxqmax(ij - iim, l) = dxqmax(ij, l) ENDDO DO ij = iip2 + 1, ip1jm IF(dxqu(ij - 1) * dxqu(ij)>0) THEN dxq(ij, l) = dxqu(ij - 1) + dxqu(ij) ELSE ! extremum local dxq(ij, l) = 0. ENDIF dxq(ij, l) = 0.5 * dxq(ij, l) dxq(ij, l) = & sign(min(abs(dxq(ij, l)), dxqmax(ij, l)), dxq(ij, l)) ENDDO ENDDO ! l=1,llm !print*,'Ok calcul des pentes' ELSE ! (pente_max.lt.-1.e-5) ! Pentes produits: ! ---------------- DO l = 1, llm DO ij = iip2, ip1jm - 1 dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq) ENDDO DO ij = iip1 + iip1, ip1jm, iip1 dxqu(ij) = dxqu(ij - iim) ENDDO DO ij = iip2 + 1, ip1jm zz(ij) = dxqu(ij - 1) * dxqu(ij) zz(ij) = zz(ij) + zz(ij) IF(zz(ij)>0) THEN dxq(ij, l) = zz(ij) / (dxqu(ij - 1) + dxqu(ij)) ELSE ! extremum local dxq(ij, l) = 0. ENDIF ENDDO ENDDO ENDIF ! (pente_max.lt.-1.e-5) ! bouclage de la pente en iip1: ! ----------------------------- DO l = 1, llm DO ij = iip1 + iip1, ip1jm, iip1 dxq(ij - iim, l) = dxq(ij, l) ENDDO DO ij = 1, ip1jmp1 iadvplus(ij, l) = 0 ENDDO ENDDO ! calcul des flux a gauche et a droite ! on cumule le flux correspondant a toutes les mailles dont la masse ! au travers de la paroi pENDant le pas de temps. DO l = 1, llm DO ij = iip2, ip1jm - 1 IF (u_m(ij, l)>0.) THEN zdum(ij, l) = 1. - u_m(ij, l) / masse(ij, l, iq) u_mq(ij, l) = u_m(ij, l) * (q(ij, l, iq) + 0.5 * zdum(ij, l) * dxq(ij, l)) ELSE zdum(ij, l) = 1. + u_m(ij, l) / masse(ij + 1, l, iq) u_mq(ij, l) = u_m(ij, l) * (q(ij + 1, l, iq) & - 0.5 * zdum(ij, l) * dxq(ij + 1, l)) ENDIF ENDDO ENDDO ! detection des points ou on advecte plus que la masse de la ! maille DO l = 1, llm DO ij = iip2, ip1jm - 1 IF(zdum(ij, l)<0) THEN iadvplus(ij, l) = 1 u_mq(ij, l) = 0. ENDIF ENDDO ENDDO DO l = 1, llm DO ij = iip1 + iip1, ip1jm, iip1 iadvplus(ij, l) = iadvplus(ij - iim, l) ENDDO ENDDO ! traitement special pour le cas ou on advecte en longitude plus que le ! contenu de la maille. ! cette partie est mal vectorisee. ! calcul du nombre de maille sur lequel on advecte plus que la maille. n0 = 0 DO l = 1, llm nl(l) = 0 DO ij = iip2, ip1jm nl(l) = nl(l) + iadvplus(ij, l) ENDDO n0 = n0 + nl(l) ENDDO IF(n0>0) THEN IF (prt_level > 2) PRINT *, & '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(iadvplus(ij, l)==1.AND.mod(ij, iip1)/=0) THEN iju = iju + 1 indu(iju) = ij ENDIF ENDDO niju = iju ! 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, iq)) u_mq(ij, l) = u_mq(ij, l) + q(ijq, l, iq) & * masse(ijq, l, iq) zu_m = zu_m - masse(ijq, l, iq) i = mod(i - 2 + iim, iim) + 1 ijq = (j - 1) * iip1 + i ENDDO ! ajout de la maille non completement advectee u_mq(ij, l) = u_mq(ij, l) + zu_m * & (q(ijq, l, iq) + 0.5 * (1. - zu_m / masse(ijq, l, iq)) & * dxq(ijq, l)) ELSE ijq = ij + 1 i = ijq - (j - 1) * iip1 ! accumulation pour les mailles completements advectees DO while(-zu_m>masse(ijq, l, iq)) u_mq(ij, l) = u_mq(ij, l) - q(ijq, l, iq) & * masse(ijq, l, iq) zu_m = zu_m + masse(ijq, l, iq) i = mod(i, iim) + 1 ijq = (j - 1) * iip1 + i ENDDO ! ajout de la maille non completement advectee u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l, iq) - & 0.5 * (1. + zu_m / masse(ijq, l, iq)) * dxq(ijq, l)) ENDIF ENDDO ENDIF ENDDO ENDIF ! n0.gt.0 ! bouclage en latitude !print*,'cvant bouclage en latitude' DO l = 1, llm DO ij = iip1 + iip1, ip1jm, iip1 u_mq(ij, l) = u_mq(ij - iim, l) ENDDO ENDDO ! CRisi: appel récursif de l'advection sur les fils. ! Il faut faire ça avant d'avoir mis à jour q et masse !WRITE(*,*) 'vlsplt 326: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen DO ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) DO l = 1, llm DO ij = iip2, ip1jm ! On a besoin de q et masse seulement entre iip2 et ip1jm !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) ! !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) !Mvals: veiller a ce qu'on n'ait pas de denominateur nul masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass) IF (q(ij, l, iq)>min_qParent) THEN Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq) else Ratio(ij, l, iq2) = min_ratio endif enddo enddo enddo DO ifils = 1, tracers(iq)%nqChildren iq2 = tracers(iq)%iqDescen(ifils) CALL vlx(Ratio, pente_max, masseq, u_mq, iq2) enddo ! end CRisi ! calcul des tENDances DO l = 1, llm DO ij = iip2 + 1, ip1jm !MVals: veiller a ce qu'on ait pas de denominateur nul new_m = max(masse(ij, l, iq) + u_m(ij - 1, l) - u_m(ij, l), min_qMass) q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + & u_mq(ij - 1, l) - u_mq(ij, l)) & / new_m masse(ij, l, iq) = new_m ENDDO DO ij = iip1 + iip1, ip1jm, iip1 q(ij - iim, l, iq) = q(ij, l, iq) masse(ij - iim, l, iq) = masse(ij, l, iq) ENDDO ENDDO ! retablir les fils en rapport de melange par rapport a l'air: ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio ! puis on boucle en longitude DO ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) DO l = 1, llm DO ij = iip2 + 1, ip1jm q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2) enddo DO ij = iip1 + iip1, ip1jm, iip1 q(ij - iim, l, iq2) = q(ij, l, iq2) enddo ! DO ij=ijb+iip1-1,ije,iip1 enddo !DO l=1,llm enddo END SUBROUTINE vlx RECURSIVE SUBROUTINE vly(q, pente_max, masse, masse_adv_v, iq) USE infotrac, ONLY: nqtot, tracers, & ! CRisi min_qParent, min_qMass, min_ratio ! MVals et CRisi USE lmdz_ssum_scopy, ONLY: ssum USE lmdz_comgeom ! Auteurs: P.Le Van, F.Hourdin, F.Forget ! ******************************************************************** ! Shema d'advection " pseudo amont " . ! ******************************************************************** ! q,masse_adv_v,w sont des arguments d'entree pour le s-pg .... ! dq sont des arguments de sortie pour le s-pg .... ! -------------------------------------------------------------------- USE comconst_mod, ONLY: pi USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE ! ! Arguments: ! ---------- REAL :: masse(ip1jmp1, llm, nqtot), pente_max REAL :: masse_adv_v(ip1jm, llm) REAL :: q(ip1jmp1, llm, nqtot) INTEGER :: iq ! CRisi ! Local ! --------- INTEGER :: i, ij, l REAL :: airej2, airejjm, airescb(iim), airesch(iim) REAL :: dyq(ip1jmp1, llm), dyqv(ip1jm) REAL :: adyqv(ip1jm), dyqmax(ip1jmp1) REAL :: qbyv(ip1jm, llm) REAL :: qpns, qpsn, dyn1, dys1, dyn2, dys2, newmasse, fn, fs LOGICAL, SAVE :: first REAL :: convpn, convps, convmpn, convmps REAL :: massepn, masseps, qpn, qps REAL :: sinlon(iip1), sinlondlon(iip1) REAL :: coslon(iip1), coslondlon(iip1) SAVE sinlon, coslon, sinlondlon, coslondlon SAVE airej2, airejjm REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) ! CRisi INTEGER :: ifils, iq2 ! CRisi DATA first/.TRUE./ ! !WRITE(*,*) 'vly 578: entree, iq=',iq IF(first) THEN PRINT*, 'Shema Amont nouveau appele dans Vanleer ' first = .FALSE. DO i = 2, iip1 coslon(i) = cos(rlonv(i)) sinlon(i) = sin(rlonv(i)) coslondlon(i) = coslon(i) * (rlonu(i) - rlonu(i - 1)) / pi sinlondlon(i) = sinlon(i) * (rlonu(i) - rlonu(i - 1)) / pi ENDDO coslon(1) = coslon(iip1) coslondlon(1) = coslondlon(iip1) sinlon(1) = sinlon(iip1) sinlondlon(1) = sinlondlon(iip1) airej2 = SSUM(iim, aire(iip2), 1) airejjm = SSUM(iim, aire(ip1jm - iim), 1) ENDIF !PRINT*,'CALCUL EN LATITUDE' DO l = 1, llm ! -------------------------------- ! CALCUL EN LATITUDE ! -------------------------------- ! On commence par calculer la valeur du traceur moyenne sur le premier cercle ! de latitude autour du pole (qpns pour le pole nord et qpsn pour ! le pole nord) qui sera utilisee pour evaluer les pentes au pole. DO i = 1, iim airescb(i) = aire(i + iip1) * q(i + iip1, l, iq) airesch(i) = aire(i + ip1jm - iip1) * q(i + ip1jm - iip1, l, iq) ENDDO qpns = SSUM(iim, airescb, 1) / airej2 qpsn = SSUM(iim, airesch, 1) / airejjm ! calcul des pentes aux points v DO ij = 1, ip1jm dyqv(ij) = q(ij, l, iq) - q(ij + iip1, l, iq) adyqv(ij) = abs(dyqv(ij)) ENDDO ! calcul des pentes aux points scalaires DO ij = iip2, ip1jm dyq(ij, l) = .5 * (dyqv(ij - iip1) + dyqv(ij)) dyqmax(ij) = min(adyqv(ij - iip1), adyqv(ij)) dyqmax(ij) = pente_max * dyqmax(ij) ENDDO ! calcul des pentes aux poles DO ij = 1, iip1 dyq(ij, l) = qpns - q(ij + iip1, l, iq) dyq(ip1jm + ij, l) = q(ip1jm + ij - iip1, l, iq) - qpsn ENDDO ! filtrage de la derivee dyn1 = 0. dys1 = 0. dyn2 = 0. dys2 = 0. DO ij = 1, iim dyn1 = dyn1 + sinlondlon(ij) * dyq(ij, l) dys1 = dys1 + sinlondlon(ij) * dyq(ip1jm + ij, l) dyn2 = dyn2 + coslondlon(ij) * dyq(ij, l) dys2 = dys2 + coslondlon(ij) * dyq(ip1jm + ij, l) ENDDO DO ij = 1, iip1 dyq(ij, l) = dyn1 * sinlon(ij) + dyn2 * coslon(ij) dyq(ip1jm + ij, l) = dys1 * sinlon(ij) + dys2 * coslon(ij) ENDDO ! calcul des pentes limites aux poles goto 8888 fn = 1. fs = 1. DO ij = 1, iim IF(pente_max * adyqv(ij)0.) THEN dyq(ij, l) = sign(min(abs(dyq(ij, l)), dyqmax(ij)), dyq(ij, l)) ELSE dyq(ij, l) = 0. ENDIF ENDDO ENDDO ! !WRITE(*,*) 'vly 756' DO l = 1, llm DO ij = 1, ip1jm IF(masse_adv_v(ij, l)>0) THEN qbyv(ij, l) = q(ij + iip1, l, iq) + dyq(ij + iip1, l) * & 0.5 * (1. - masse_adv_v(ij, l) & / masse(ij + iip1, l, iq)) ELSE qbyv(ij, l) = q(ij, l, iq) - dyq(ij, l) * & 0.5 * (1. + masse_adv_v(ij, l) & / masse(ij, l, iq)) ENDIF qbyv(ij, l) = masse_adv_v(ij, l) * qbyv(ij, l) ENDDO ENDDO ! CRisi: appel récursif de l'advection sur les fils. ! Il faut faire ça avant d'avoir mis à jour q et masse ! WRITE(*,*) 'vly 689: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen DO ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) DO l = 1, llm DO ij = 1, ip1jmp1 ! ! attention, chaque fils doit avoir son masseq, sinon, le 1er ! ! fils ecrase le masseq de ses freres. ! !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) ! !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass) IF (q(ij, l, iq)>min_qParent) THEN Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq) else Ratio(ij, l, iq2) = min_ratio endif enddo enddo enddo DO ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) CALL vly(Ratio, pente_max, masseq, qbyv, iq2) enddo DO l = 1, llm DO ij = iip2, ip1jm newmasse = masse(ij, l, iq) & + masse_adv_v(ij, l) - masse_adv_v(ij - iip1, l) q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + qbyv(ij, l) & - qbyv(ij - iip1, l)) / newmasse masse(ij, l, iq) = newmasse ENDDO convpn = SSUM(iim, qbyv(1, l), 1) convmpn = ssum(iim, masse_adv_v(1, l), 1) massepn = ssum(iim, masse(1, l, iq), 1) qpn = 0. DO ij = 1, iim qpn = qpn + masse(ij, l, iq) * q(ij, l, iq) enddo qpn = (qpn + convpn) / (massepn + convmpn) DO ij = 1, iip1 q(ij, l, iq) = qpn enddo convps = -SSUM(iim, qbyv(ip1jm - iim, l), 1) convmps = -ssum(iim, masse_adv_v(ip1jm - iim, l), 1) masseps = ssum(iim, masse(ip1jm + 1, l, iq), 1) qps = 0. DO ij = ip1jm + 1, ip1jmp1 - 1 qps = qps + masse(ij, l, iq) * q(ij, l, iq) enddo qps = (qps + convps) / (masseps + convmps) DO ij = ip1jm + 1, ip1jmp1 q(ij, l, iq) = qps enddo ENDDO ! retablir les fils en rapport de melange par rapport a l'air: DO ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) DO l = 1, llm DO ij = 1, ip1jmp1 q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2) enddo enddo enddo ! !WRITE(*,*) 'vly 853: sortie' END SUBROUTINE vly RECURSIVE SUBROUTINE vlz(q, pente_max, masse, w, iq) USE infotrac, ONLY: nqtot, tracers, & ! CRisi min_qParent, min_qMass, min_ratio ! MVals et CRisi ! Auteurs: P.Le Van, F.Hourdin, F.Forget ! ******************************************************************** ! Shema d'advection " pseudo amont " . ! ******************************************************************** ! q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... ! dq sont des arguments de sortie pour le s-pg .... ! -------------------------------------------------------------------- USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE ! ! Arguments: ! ---------- REAL :: masse(ip1jmp1, llm, nqtot), pente_max REAL :: q(ip1jmp1, llm, nqtot) REAL :: w(ip1jmp1, llm + 1) INTEGER :: iq ! Local ! --------- INTEGER :: ij, l REAL :: wq(ip1jmp1, llm + 1), newmasse REAL :: dzq(ip1jmp1, llm), dzqw(ip1jmp1, llm), adzqw(ip1jmp1, llm), dzqmax REAL :: sigw REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) ! CRisi INTEGER :: ifils, iq2 ! CRisi #ifdef BIDON REAL :: temps0,temps1,second SAVE temps0,temps1 DATA temps0,temps1/0.,0./ #endif ! On oriente tout dans le sens de la pression c'est a dire dans le ! sens de W DO l = 2, llm DO ij = 1, ip1jmp1 dzqw(ij, l) = q(ij, l - 1, iq) - q(ij, l, iq) adzqw(ij, l) = abs(dzqw(ij, l)) ENDDO ENDDO DO l = 2, llm - 1 DO ij = 1, ip1jmp1 IF(dzqw(ij, l) * dzqw(ij, l + 1)>0.) THEN dzq(ij, l) = 0.5 * (dzqw(ij, l) + dzqw(ij, l + 1)) ELSE dzq(ij, l) = 0. ENDIF dzqmax = pente_max * min(adzqw(ij, l), adzqw(ij, l + 1)) dzq(ij, l) = sign(min(abs(dzq(ij, l)), dzqmax), dzq(ij, l)) ENDDO ENDDO ! !WRITE(*,*) 'vlz 954' DO ij = 1, ip1jmp1 dzq(ij, 1) = 0. dzq(ij, llm) = 0. ENDDO ! --------------------------------------------------------------- ! .... calcul des termes d'advection verticale ....... ! --------------------------------------------------------------- ! calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq DO l = 1, llm - 1 DO ij = 1, ip1jmp1 IF(w(ij, l + 1)>0.) THEN sigw = w(ij, l + 1) / masse(ij, l + 1, iq) wq(ij, l + 1) = w(ij, l + 1) * (q(ij, l + 1, iq) & + 0.5 * (1. - sigw) * dzq(ij, l + 1)) ELSE sigw = w(ij, l + 1) / masse(ij, l, iq) wq(ij, l + 1) = w(ij, l + 1) * (q(ij, l, iq) - 0.5 * (1. + sigw) * dzq(ij, l)) ENDIF ENDDO ENDDO DO ij = 1, ip1jmp1 wq(ij, llm + 1) = 0. wq(ij, 1) = 0. ENDDO ! CRisi: appel récursif de l'advection sur les fils. ! Il faut faire ça avant d'avoir mis à jour q et masse ! !WRITE(*,*) 'vlsplt 942: iq,nqChildren(iq)=',iq,nqChildren(iq) DO ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) DO l = 1, llm DO ij = 1, ip1jmp1 ! !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) ! !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass) IF (q(ij, l, iq)>min_qParent) THEN Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq) else Ratio(ij, l, iq2) = min_ratio endif enddo enddo enddo DO ifils = 1, tracers(iq)%nqChildren iq2 = tracers(iq)%iqDescen(ifils) CALL vlz(Ratio, pente_max, masseq, wq, iq2) enddo ! end CRisi DO l = 1, llm DO ij = 1, ip1jmp1 newmasse = masse(ij, l, iq) + w(ij, l + 1) - w(ij, l) q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + wq(ij, l + 1) - wq(ij, l)) & / newmasse masse(ij, l, iq) = newmasse ENDDO ENDDO ! retablir les fils en rapport de melange par rapport a l'air: DO ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) DO l = 1, llm DO ij = 1, ip1jmp1 q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2) enddo enddo enddo END SUBROUTINE vlz SUBROUTINE minmaxq(zq, qmin, qmax, comment) USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet CHARACTER(LEN = 20) :: comment REAL :: qmin, qmax REAL :: zq(ip1jmp1, llm) REAL :: zzq(iip1, jjp1, llm) END SUBROUTINE minmaxq