! $Id: vlsplt_loc.f90 5128 2024-07-25 15:47:25Z abarral $ SUBROUTINE vlx_loc(q, pente_max, masse, u_m, ijb_x, ije_x, iq) ! 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 parallel_lmdz USE infotrac, ONLY: nqtot, tracers, & ! CRisi & min_qParent, min_qMass, min_ratio ! MVals et CRisi USE lmdz_iniprint, ONLY: lunout, prt_level IMPLICIT NONE ! include "dimensions.h" include "paramet.h" ! ! ! Arguments: ! ---------- REAL :: masse(ijb_u:ije_u, llm, nqtot), pente_max REAL :: u_m(ijb_u:ije_u, llm), pbarv(iip1, jjb_v:jje_v, llm) REAL :: q(ijb_u:ije_u, llm, nqtot) ! CRisi: ajout dimension nqtot REAL :: w(ijb_u:ije_u, llm) INTEGER :: iq ! CRisi ! ! Local ! --------- ! INTEGER :: ij, l, j, i, iju, ijq, indu(ijnb_u), niju INTEGER :: n0, iadvplus(ijb_u:ije_u, llm), nl(llm) ! REAL :: new_m, zu_m, zdum(ijb_u:ije_u, llm) REAL :: sigu(ijb_u:ije_u), dxq(ijb_u:ije_u, llm), dxqu(ijb_u:ije_u) REAL :: zz(ijb_u:ije_u) REAL :: adxqu(ijb_u:ije_u), dxqmax(ijb_u:ije_u, llm) REAL :: u_mq(ijb_u:ije_u, llm) REAL :: Ratio(ijb_u:ije_u, llm, nqtot) ! CRisi INTEGER :: ifils, iq2 ! CRisi Logical :: extremum REAL :: z1, z2, z3 INTEGER :: ijb, ije, ijb_x, ije_x !WRITE(*,*) 'vlsplt 58: entree dans vlx_loc, iq,ijb_x=', ! & iq,ijb_x ! calcul de la pente a droite et a gauche de la maille ijb = ijb_x ije = ije_x IF (pole_nord.AND.ijb==1) ijb = ijb + iip1 IF (pole_sud.AND.ije==ip1jmp1) ije = ije - iip1 IF (pente_max>-1.e-5) THEN ! IF (pente_max.gt.10) THEN ! calcul des pentes avec limitation, Van Leer scheme I: ! ----------------------------------------------------- ! on a besoin de q entre ijb et ije ! calcul de la pente aux points u !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb, ije - 1 dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq) ! IF(u_m(ij,l).lt.0.) stop 'limx n admet pas les U<0' ! sigu(ij)=u_m(ij,l)/masse(ij,l,iq) ENDDO DO ij = ijb + iip1 - 1, ije, iip1 dxqu(ij) = dxqu(ij - iim) ! sigu(ij)=sigu(ij-iim) ENDDO DO ij = ijb, ije adxqu(ij) = abs(dxqu(ij)) ENDDO ! calcul de la pente maximum dans la maille en valeur absolue DO ij = ijb + 1, ije 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 = ijb + iip1 - 1, ije, iip1 dxqmax(ij - iim, l) = dxqmax(ij, l) ENDDO DO ij = ijb + 1, ije 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 !$OMP END DO NOWAIT ! PRINT*,'Ok calcul des pentes' ELSE ! (pente_max.lt.-1.e-5) ! Pentes produits: ! ---------------- !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb, ije - 1 dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq) ENDDO DO ij = ijb + iip1 - 1, ije, iip1 dxqu(ij) = dxqu(ij - iim) ENDDO DO ij = ijb + 1, ije 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 !$OMP END DO NOWAIT ENDIF ! (pente_max.lt.-1.e-5) !WRITE(*,*) 'vlx 156: iq,ijb_x=',iq,ijb_x ! bouclage de la pente en iip1: ! ----------------------------- !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb + iip1 - 1, ije, iip1 dxq(ij - iim, l) = dxq(ij, l) ENDDO DO ij = ijb, ije iadvplus(ij, l) = 0 ENDDO ENDDO !$OMP END DO NOWAIT ! PRINT*,'Bouclage en iip1' ! 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. ! PRINT*,'Cumule ....' !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) ! on a besoin de masse entre ijb et ije DO l = 1, llm DO ij = ijb, ije - 1 ! PRINT*,'masse(',ij,')=',masse(ij,l,iq) 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 !$OMP END DO NOWAIT ! detection des points ou on advecte plus que la masse de la ! maille !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb, ije - 1 IF(zdum(ij, l)<0) THEN iadvplus(ij, l) = 1 u_mq(ij, l) = 0. ENDIF ENDDO ENDDO !$OMP END DO NOWAIT ! PRINT*,'Ok test 1' !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb + iip1 - 1, ije, iip1 iadvplus(ij, l) = iadvplus(ij - iim, l) ENDDO ENDDO !$OMP END DO NOWAIT ! PRINT*,'Ok test 2' ! 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 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm nl(l) = 0 DO ij = ijb, ije nl(l) = nl(l) + iadvplus(ij, l) ENDDO n0 = n0 + nl(l) ENDDO !$OMP END DO NOWAIT !ym IF(n0.gt.1) THEN !ym IF(n0.gt.0) THEN ! PRINT*,'Nombre de points pour lesquels on advect plus que le' ! & ,'contenu de la maille : ',n0 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm IF(nl(l)>0) THEN iju = 0 ! indicage des mailles concernees par le traitement special DO ij = ijb, ije IF(iadvplus(ij, l)==1.AND.mod(ij, iip1)/=0) THEN iju = iju + 1 indu(iju) = ij ENDIF ENDDO niju = iju !PRINT*,'vlx 278, 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, 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 !$OMP END DO NOWAIT !ym ENDIF ! n0.gt.0 ! bouclage en latitude ! PRINT*,'Avant bouclage en latitude' !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb + iip1 - 1, ije, iip1 u_mq(ij, l) = u_mq(ij - iim, l) ENDDO ENDDO !$OMP END DO NOWAIT ! CRisi: appel récursif de l'advection sur les fils. ! Il faut faire ça avant d'avoir mis à jour q et masse do ifils = 1, tracers(iq)%nqDescen ! attention: comme Ratio est utilisé comme q dans l'appel ! recursif, il doit contenir à lui seul tous les indices de tous ! les descendants! iq2 = tracers(iq)%iqDescen(ifils) !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb, ije ! On a besoin de q et masse seulement entre ijb et ije. On ne ! les calcule donc que de ijb à ije !MVals: veiller a ce qu'on n'ait pas de denominateur nul masse(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass) IF (q(ij, l, iq)>min_qParent) then ! modif 13 nov 2020 Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq) else Ratio(ij, l, iq2) = min_ratio endif enddo enddo !$OMP END DO NOWAIT enddo !do ifils=1,tracers(iq)%nqDescen do ifils = 1, tracers(iq)%nqChildren iq2 = tracers(iq)%iqDescen(ifils) CALL vlx_loc(Ratio, pente_max, masse, u_mq, ijb_x, ije_x, iq2) enddo ! end CRisi ! calcul des tENDances !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb + 1, ije !MVals: veiller a ce qu'on n'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 ! ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous) DO ij = ijb + iip1 - 1, ije, iip1 q(ij - iim, l, iq) = q(ij, l, iq) masse(ij - iim, l, iq) = masse(ij, l, iq) ENDDO ENDDO !$OMP END DO NOWAIT ! retablir les fils en rapport de melange par rapport a l'air: ! On calcule q entre ijb+1 et ije -> on fait pareil pour ratio ! puis on boucle en longitude do ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb + 1, ije q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2) enddo DO ij = ijb + iip1 - 1, ije, iip1 q(ij - iim, l, iq2) = q(ij, l, iq2) enddo enddo !$OMP END DO NOWAIT enddo END SUBROUTINE vlx_loc SUBROUTINE vly_loc(q, pente_max, masse, masse_adv_v, iq) ! ! 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 parallel_lmdz USE infotrac, ONLY: nqtot, tracers, & ! CRisi & min_qParent, min_qMass, min_ratio ! MVals et CRisi USE comconst_mod, ONLY: pi USE lmdz_ssum_scopy, ONLY: ssum IMPLICIT NONE ! include "dimensions.h" include "paramet.h" include "comgeom.h" ! ! ! Arguments: ! ---------- REAL :: masse(ijb_u:ije_u, llm, nqtot), pente_max REAL :: masse_adv_v(ijb_v:ije_v, llm) REAL :: q(ijb_u:ije_u, llm, nqtot), dq(ijb_u:ije_u, llm) INTEGER :: iq ! CRisi ! ! Local ! --------- ! INTEGER :: i, ij, l ! REAL :: airej2, airejjm, airescb(iim), airesch(iim) REAL :: dyq(ijb_u:ije_u, llm), dyqv(ijb_v:ije_v), zdvm(ijb_u:ije_u, llm) REAL :: adyqv(ijb_v:ije_v), dyqmax(ijb_u:ije_u) REAL :: qbyv(ijb_v:ije_v, llm) REAL :: qpns, qpsn, appn, apps, dyn1, dys1, dyn2, dys2, newmasse, fn, fs ! REAL newq,oldmasse Logical :: extremum, first REAL :: temps0, temps1, temps2, temps3, temps4, temps5, second SAVE temps0, temps1, temps2, temps3, temps4, temps5 !$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5) SAVE first !$OMP THREADPRIVATE(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 !$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon) SAVE airej2, airejjm !$OMP THREADPRIVATE(airej2,airejjm) REAL :: Ratio(ijb_u:ije_u, llm, nqtot) ! CRisi INTEGER :: ifils, iq2 ! CRisi DATA first/.TRUE./ DATA temps0, temps1, temps2, temps3, temps4, temps5/0., 0., 0., 0., 0., 0./ INTEGER :: ijb, ije INTEGER :: ijbm, ijem ijb = ij_begin - 2 * iip1 ije = ij_end + 2 * iip1 IF (pole_nord) ijb = ij_begin IF (pole_sud) ije = ij_end 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' !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 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. IF (pole_nord) THEN DO i = 1, iim airescb(i) = aire(i + iip1) * q(i + iip1, l, iq) ENDDO qpns = SSUM(iim, airescb, 1) / airej2 endif IF (pole_sud) THEN DO i = 1, iim airesch(i) = aire(i + ip1jm - iip1) * q(i + ip1jm - iip1, l, iq) ENDDO qpsn = SSUM(iim, airesch, 1) / airejjm endif ! calcul des pentes aux points v ijb = ij_begin - 2 * iip1 ije = ij_end + iip1 IF (pole_nord) ijb = ij_begin IF (pole_sud) ije = ij_end - iip1 ! on a besoin de q entre ij_begin-2*iip1 et ij_end+2*iip1 ! Si pole sud, entre ij_begin-2*iip1 et ij_end ! Si pole Nord, entre ij_begin et ij_end+2*iip1 DO ij = ijb, ije dyqv(ij) = q(ij, l, iq) - q(ij + iip1, l, iq) adyqv(ij) = abs(dyqv(ij)) ENDDO ! calcul des pentes aux points scalaires ijb = ij_begin - iip1 ije = ij_end + iip1 IF (pole_nord) ijb = ij_begin + iip1 IF (pole_sud) ije = ij_end - iip1 DO ij = ijb, ije 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 IF (pole_nord) THEN DO ij = 1, iip1 dyq(ij, l) = qpns - q(ij + iip1, l, iq) ENDDO dyn1 = 0. dyn2 = 0. DO ij = 1, iim dyn1 = dyn1 + sinlondlon(ij) * dyq(ij, l) dyn2 = dyn2 + coslondlon(ij) * dyq(ij, l) ENDDO DO ij = 1, iip1 dyq(ij, l) = dyn1 * sinlon(ij) + dyn2 * coslon(ij) ENDDO DO ij = 1, iip1 dyq(ij, l) = 0. ENDDO ! ym tout cela ne sert pas a grand chose ENDIF IF (pole_sud) THEN DO ij = 1, iip1 dyq(ip1jm + ij, l) = q(ip1jm + ij - iip1, l, iq) - qpsn ENDDO dys1 = 0. dys2 = 0. DO ij = 1, iim dys1 = dys1 + sinlondlon(ij) * dyq(ip1jm + ij, l) dys2 = dys2 + coslondlon(ij) * dyq(ip1jm + ij, l) ENDDO DO ij = 1, iip1 dyq(ip1jm + ij, l) = dys1 * sinlon(ij) + dys2 * coslon(ij) ENDDO DO ij = 1, iip1 dyq(ip1jm + ij, l) = 0. ENDDO ! ym tout cela ne sert pas a grand chose ENDIF ! filtrage de la derivee ! calcul des pentes limites aux poles ! ym partie inutile ! goto 8888 ! fn=1. ! fs=1. ! DO ij=1,iim ! IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN ! fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn) ! ENDIF ! IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN ! fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs) ! ENDIF ! ENDDO ! DO ij=1,iip1 ! dyq(ij,l)=fn*dyq(ij,l) ! dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l) ! ENDDO ! 8888 continue !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! En memoire de dIFferents tests sur la ! limitation des pentes aux poles. !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! PRINT*,dyq(1) ! PRINT*,dyqv(iip1+1) ! appn=abs(dyq(1)/dyqv(iip1+1)) ! PRINT*,dyq(ip1jm+1) ! PRINT*,dyqv(ip1jm-iip1+1) ! apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) ! DO ij=2,iim ! appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) ! apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) ! ENDDO ! appn=min(pente_max/appn,1.) ! apps=min(pente_max/apps,1.) ! ! ! cas ou on a un extremum au pole ! ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) ! & appn=0. ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* ! & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) ! & apps=0. ! ! limitation des pentes aux poles ! DO ij=1,iip1 ! dyq(ij)=appn*dyq(ij) ! dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) ! ENDDO ! ! test ! DO ij=1,iip1 ! dyq(iip1+ij)=0. ! dyq(ip1jm+ij-iip1)=0. ! ENDDO ! DO ij=1,ip1jmp1 ! dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1)) ! ENDDO ! ! changement 10 07 96 ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) ! & THEN ! DO ij=1,iip1 ! dyqmax(ij)=0. ! ENDDO ! ELSE ! DO ij=1,iip1 ! dyqmax(ij)=pente_max*abs(dyqv(ij)) ! ENDDO ! ENDIF ! ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* ! & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) ! &THEN ! DO ij=ip1jm+1,ip1jmp1 ! dyqmax(ij)=0. ! ENDDO ! ELSE ! DO ij=ip1jm+1,ip1jmp1 ! dyqmax(ij)=pente_max*abs(dyqv(ij-iip1)) ! ENDDO ! ENDIF ! fin changement 10 07 96 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! calcul des pentes limitees ijb = ij_begin - iip1 ije = ij_end + iip1 IF (pole_nord) ijb = ij_begin + iip1 IF (pole_sud) ije = ij_end - iip1 DO ij = ijb, ije IF(dyqv(ij) * dyqv(ij - iip1)>0.) THEN dyq(ij, l) = sign(min(abs(dyq(ij, l)), dyqmax(ij)), dyq(ij, l)) ELSE dyq(ij, l) = 0. ENDIF ENDDO ENDDO !$OMP END DO NOWAIT ijb = ij_begin - iip1 ije = ij_end IF (pole_nord) ijb = ij_begin IF (pole_sud) ije = ij_end - iip1 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb, ije 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 !$OMP END DO NOWAIT ! 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,nqChildren(iq)=',iq,tracers(iq)%nqChildren ijb = ij_begin - 2 * iip1 ije = ij_end + 2 * iip1 ijbm = ij_begin - iip1 ijem = ij_end + iip1 IF (pole_nord) ijb = ij_begin IF (pole_sud) ije = ij_end IF (pole_nord) ijbm = ij_begin IF (pole_sud) ijem = ij_end do ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm ! modif des bornes: CRisi 16 nov 2020 ! d'abord masse avec bornes corrigées DO ij = ijbm, ijem !MVals: veiller a ce qu'on n'ait pas de denominateur nul masse(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass) enddo ! ensuite Ratio avec anciennes bornes DO ij = ijb, ije !MVals: veiller a ce qu'on n'ait pas de denominateur nul IF (q(ij, l, iq)>min_qParent) then ! modif 13 nov 2020 Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq) else Ratio(ij, l, iq2) = min_ratio endif enddo !DO ij=ijbm,ijem enddo !DO l=1,llm !$OMP END DO NOWAIT enddo do ifils = 1, tracers(iq)%nqChildren iq2 = tracers(iq)%iqDescen(ifils) CALL vly_loc(Ratio, pente_max, masse, qbyv, iq2) enddo ! end CRisi ijb = ij_begin ije = ij_end IF (pole_nord) ijb = ij_begin + iip1 IF (pole_sud) ije = ij_end - iip1 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb, ije 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 IF (pole_nord) THEN 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 endif IF (pole_sud) THEN 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 endif ENDDO !$OMP END DO NOWAIT ! retablir les fils en rapport de melange par rapport a l'air: ijb = ij_begin ije = ij_end ! if (pole_nord) ijb=ij_begin ! if (pole_sud) ije=ij_end do ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb, ije q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2) enddo enddo !$OMP END DO NOWAIT enddo END SUBROUTINE vly_loc SUBROUTINE vlz_loc(q, pente_max, masse, w, ijb_x, ije_x, iq) ! ! 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 parallel_lmdz USE vlz_mod USE infotrac, ONLY: nqtot, tracers, & ! CRisi & min_qParent, min_qMass, min_ratio ! MVals et CRisi USE lmdz_iniprint, ONLY: lunout, prt_level IMPLICIT NONE ! include "dimensions.h" include "paramet.h" ! ! ! Arguments: ! ---------- REAL :: masse(ijb_u:ije_u, llm, nqtot), pente_max REAL :: q(ijb_u:ije_u, llm, nqtot) REAL :: w(ijb_u:ije_u, llm + 1, nqtot) INTEGER :: iq ! ! Local ! --------- ! INTEGER :: i, ij, l, j, ii REAL, DIMENSION(ijb_u:ije_u, llm + 1) :: wresi, morig, qorig, dzqorig INTEGER, DIMENSION(ijb_u:ije_u, llm + 1) :: lorig INTEGER, SAVE :: countcfl !$OMP THREADPRIVATE(countcfl) ! REAL :: newmasse REAL :: dzqmax REAL :: sigw REAL :: temps0, temps1, temps2, temps3, temps4, temps5, second SAVE temps0, temps1, temps2, temps3, temps4, temps5 !$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5) DATA temps0, temps1, temps2, temps3, temps4, temps5/0., 0., 0., 0., 0., 0./ INTEGER :: ijb, ije, ijb_x, ije_x LOGICAL, SAVE :: first = .TRUE. !$OMP THREADPRIVATE(first) !REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi ! Ces varibles doivent être déclarées en pointer et en save dans ! vlz_loc si on veut qu'elles soient vues par tous les threads. INTEGER :: ifils, iq2 ! CRisi IF (first) THEN first = .FALSE. ENDIF ! On oriente tout dans le sens de la pression c'est a dire dans le ! sens de W !WRITE(*,*) 'vlsplt 926: entree dans vlz_loc, iq=',iq ijb = ijb_x ije = ije_x !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 2, llm DO ij = ijb, ije dzqw(ij, l) = q(ij, l - 1, iq) - q(ij, l, iq) adzqw(ij, l) = abs(dzqw(ij, l)) ENDDO ENDDO !$OMP END DO !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 2, llm - 1 DO ij = ijb, ije 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 !$OMP END DO NOWAIT !$OMP MASTER DO ij = ijb, ije dzq(ij, 1) = 0. dzq(ij, llm) = 0. ENDDO !$OMP END MASTER !$OMP BARRIER !-------------------------------------------------------- ! On repere les points qui violent le CFL (|w| > masse) !-------------------------------------------------------- countcfl = 0 ! PRINT*,'vlz nouveau' !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 2, llm DO ij = ijb, ije IF((w(ij, l, iq)>0.AND.w(ij, l, iq)>masse(ij, l, iq)) & .OR. (w(ij, l, iq)<=0.AND.ABS(w(ij, l, iq))>masse(ij, l - 1, iq))) & countcfl = countcfl + 1 ENDDO ENDDO !$OMP END DO NOWAIT ! --------------------------------------------------------------- ! Identification des mailles ou on viole le CFL : w > masse ! --------------------------------------------------------------- IF (countcfl==0) THEN ! --------------------------------------------------------------- ! .... calcul des termes d'advection verticale ....... ! Dans le cas où le |w| < masse partout. ! Version d'origine ! Pourrait etre enleve si on voit que le code plus general ! est aussi rapide ! --------------------------------------------------------------- ! calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq ! !WRITE(*,*) 'vlz 982,ijb,ije=',ijb,ije !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm - 1 do ij = ijb, ije IF(w(ij, l + 1, iq)>0.) THEN sigw = w(ij, l + 1, iq) / masse(ij, l + 1, iq) wq(ij, l + 1, iq) = w(ij, l + 1, iq) * (q(ij, l + 1, iq) & + 0.5 * (1. - sigw) * dzq(ij, l + 1)) ELSE sigw = w(ij, l + 1, iq) / masse(ij, l, iq) wq(ij, l + 1, iq) = w(ij, l + 1, iq) * (q(ij, l, iq) & - 0.5 * (1. + sigw) * dzq(ij, l)) ENDIF ENDDO ENDDO !$OMP END DO NOWAIT !WRITE(*,*) 'vlz 1001' ELSE ! countcfl>=1 IF (prt_level>9) THEN WRITE(lunout, *)'vlz passage dans le non local' ENDIF ! --------------------------------------------------------------- ! Debut du traitement du cas ou on viole le CFL : w > masse ! --------------------------------------------------------------- ! Initialisation !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 2, llm DO ij = ijb, ije wresi(ij, l) = w(ij, l, iq) wq(ij, l, iq) = 0. IF(w(ij, l, iq)>0.) THEN lorig(ij, l) = l morig(ij, l) = masse(ij, l, iq) qorig(ij, l) = q(ij, l, iq) dzqorig(ij, l) = dzq(ij, l) ELSE lorig(ij, l) = l - 1 morig(ij, l) = masse(ij, l - 1, iq) qorig(ij, l) = q(ij, l - 1, iq) dzqorig(ij, l) = dzq(ij, l - 1) ENDIF ENDDO ENDDO !$OMP END DO NOWAIT ! Reindicage vertical en accumulant les flux sur ! les mailles qui viollent le CFL ! on itère jusqu'à ce que tous les poins satisfassent ! le critère DO WHILE (countcfl>=1) IF (prt_level>9) THEN WRITE(lunout, *)'On viole le CFL Vertical sur ', countcfl, ' pts' ENDIF countcfl = 0 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 2, llm DO ij = ijb, ije IF (ABS(wresi(ij, l))>morig(ij, l)) THEN countcfl = countcfl + 1 ! rm : les 8 lignes ci dessous pourraient sans doute s'ecrire ! avec la fonction sign IF(w(ij, l, iq)>0.) THEN wresi(ij, l) = wresi(ij, l) - morig(ij, l) wq(ij, l, iq) = wq(ij, l, iq) + morig(ij, l) * qorig(ij, l) lorig(ij, l) = lorig(ij, l) + 1 ELSE wresi(ij, l) = wresi(ij, l) + morig(ij, l) wq(ij, l, iq) = wq(ij, l, iq) - morig(ij, l) * qorig(ij, l) lorig(ij, l) = lorig(ij, l) - 1 ENDIF ! CRisi 24nov2020: ajout d'un message d'erreur clair au lieu d'un plantage ! pour seg fault IF (lorig(ij, l)==0) THEN CALL abort_gcm("vlz in vlsplt_loc", & "unfixable violation of CFL", 1) endif morig(ij, l) = masse(ij, lorig(ij, l), iq) qorig(ij, l) = q(ij, lorig(ij, l), iq) dzqorig(ij, l) = dzq(ij, lorig(ij, l)) ENDIF ENDDO ENDDO !$OMP END DO NOWAIT ENDDO ! WHILE (countcfl>=1) !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 2, llm do ij = ijb, ije sigw = wresi(ij, l) / morig(ij, l) IF(w(ij, l, iq)>0.) THEN wq(ij, l, iq) = wq(ij, l, iq) + wresi(ij, l) * (qorig(ij, l) & + 0.5 * (1. - sigw) * dzqorig(ij, l)) ELSE wq(ij, l, iq) = wq(ij, l, iq) + wresi(ij, l) * (qorig(ij, l) & - 0.5 * (1. + sigw) * dzqorig(ij, l)) ENDIF ENDDO ENDDO !$OMP END DO NOWAIT ENDIF ! councfl=0 !$OMP MASTER DO ij = ijb, ije wq(ij, llm + 1, iq) = 0. wq(ij, 1, iq) = 0. ENDDO !$OMP END MASTER !$OMP BARRIER ! 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,tracers(iq)%nqChildren do ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb, ije !MVals: veiller a ce qu'on n'ait pas de denominateur nul masse(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 !wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015 w(ij, l, iq2) = wq(ij, l, iq) enddo enddo !$OMP END DO NOWAIT enddo !$OMP BARRIER do ifils = 1, tracers(iq)%nqChildren iq2 = tracers(iq)%iqDescen(ifils) CALL vlz_loc(Ratio, pente_max, masse, w, ijb_x, ije_x, iq2) enddo ! end CRisi ! CRisi: On rajoute ici une barrière car on veut être sur que tous les ! wq soient synchronisés !$OMP BARRIER !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb, ije newmasse = masse(ij, l, iq) + w(ij, l + 1, iq) - w(ij, l, iq) q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) & + wq(ij, l + 1, iq) - wq(ij, l, iq)) & / newmasse masse(ij, l, iq) = newmasse ENDDO ENDDO !$OMP END DO NOWAIT ! retablir les fils en rapport de melange par rapport a l'air: do ifils = 1, tracers(iq)%nqDescen iq2 = tracers(iq)%iqDescen(ifils) !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb, ije q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2) enddo enddo !$OMP END DO NOWAIT enddo END SUBROUTINE vlz_loc