Changeset 5248 for LMDZ6/trunk/libf/dyn3dmem/vlsplt_loc.F90
- Timestamp:
- Oct 21, 2024, 7:05:31 PM (17 hours ago)
- File:
-
- 1 moved
-
LMDZ6/trunk/libf/dyn3dmem/vlsplt_loc.F90 (moved) (moved from LMDZ6/trunk/libf/dyn3dmem/vlsplt_loc.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/dyn3dmem/vlsplt_loc.F90
r5247 r5248 2 2 ! $Id$ 3 3 ! 4 5 6 c Auteurs: P.Le Van, F.Hourdin, F.Forget 7 c 8 c********************************************************************9 cShema d'advection " pseudo amont " .10 c********************************************************************11 cnq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....12 c 13 c 14 c--------------------------------------------------------------------15 16 USE infotrac, ONLY : nqtot,tracers,! CRisi &17 &min_qParent,min_qMass,min_ratio ! MVals et CRisi18 19 c 20 21 22 23 c 24 c 25 cArguments:26 c----------27 REALmasse(ijb_u:ije_u,llm,nqtot),pente_max28 REALu_m( ijb_u:ije_u,llm),pbarv( iip1,jjb_v:jje_v,llm)29 REALq(ijb_u:ije_u,llm,nqtot) ! CRisi: ajout dimension nqtot30 REAL w(ijb_u:ije_u,llm)31 INTEGERiq ! CRisi32 c 33 c Local 34 c---------35 c 36 INTEGERij,l,j,i,iju,ijq,indu(ijnb_u),niju37 INTEGERn0,iadvplus(ijb_u:ije_u,llm),nl(llm)38 c 39 REALnew_m,zu_m,zdum(ijb_u:ije_u,llm)40 REALsigu(ijb_u:ije_u),dxq(ijb_u:ije_u,llm),dxqu(ijb_u:ije_u)41 REALzz(ijb_u:ije_u)42 REALadxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm)43 REALu_mq(ijb_u:ije_u,llm)44 45 REALRatio(ijb_u:ije_u,llm,nqtot) ! CRisi46 INTEGERifils,iq2 ! CRisi47 48 Logicalextremum49 50 REALSSUM51 52 53 REALz1,z2,z354 55 INTEGERijb,ije,ijb_x,ije_x56 57 58 !& iq,ijb_x59 ccalcul de la pente a droite et a gauche de la maille60 61 62 63 64 65 66 67 68 cIF (pente_max.gt.10) THEN69 70 ccalcul des pentes avec limitation, Van Leer scheme I:71 c-----------------------------------------------------72 73 ccalcul de la pente aux points u74 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 75 76 77 78 79 cIF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'80 csigu(ij)=u_m(ij,l)/masse(ij,l,iq)81 82 83 84 csigu(ij)=sigu(ij-iim)85 86 87 88 89 90 91 ccalcul de la pente maximum dans la maille en valeur absolue92 93 94 dxqmax(ij,l)=pente_max*95 ,min(adxqu(ij-1),adxqu(ij))96 climitation subtile97 c, min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))98 99 100 101 102 103 104 105 106 4 RECURSIVE SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x,iq) 5 6 ! Auteurs: P.Le Van, F.Hourdin, F.Forget 7 ! 8 ! ******************************************************************** 9 ! Shema d'advection " pseudo amont " . 10 ! ******************************************************************** 11 ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... 12 ! 13 ! 14 ! -------------------------------------------------------------------- 15 USE parallel_lmdz 16 USE infotrac, ONLY : nqtot,tracers, & ! CRisi & 17 min_qParent,min_qMass,min_ratio ! MVals et CRisi 18 IMPLICIT NONE 19 ! 20 include "dimensions.h" 21 include "paramet.h" 22 include "iniprint.h" 23 ! 24 ! 25 ! Arguments: 26 ! ---------- 27 REAL :: masse(ijb_u:ije_u,llm,nqtot),pente_max 28 REAL :: u_m( ijb_u:ije_u,llm),pbarv( iip1,jjb_v:jje_v,llm) 29 REAL :: q(ijb_u:ije_u,llm,nqtot) ! CRisi: ajout dimension nqtot 30 REAL :: w(ijb_u:ije_u,llm) 31 INTEGER :: iq ! CRisi 32 ! 33 ! Local 34 ! --------- 35 ! 36 INTEGER :: ij,l,j,i,iju,ijq,indu(ijnb_u),niju 37 INTEGER :: n0,iadvplus(ijb_u:ije_u,llm),nl(llm) 38 ! 39 REAL :: new_m,zu_m,zdum(ijb_u:ije_u,llm) 40 REAL :: sigu(ijb_u:ije_u),dxq(ijb_u:ije_u,llm),dxqu(ijb_u:ije_u) 41 REAL :: zz(ijb_u:ije_u) 42 REAL :: adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm) 43 REAL :: u_mq(ijb_u:ije_u,llm) 44 45 REAL :: Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi 46 INTEGER :: ifils,iq2 ! CRisi 47 48 Logical :: extremum 49 50 REAL :: SSUM 51 EXTERNAL SSUM 52 53 REAL :: z1,z2,z3 54 55 INTEGER :: ijb,ije,ijb_x,ije_x 56 57 ! !write(*,*) 'vlsplt 58: entree dans vlx_loc, iq,ijb_x=', 58 ! & iq,ijb_x 59 ! calcul de la pente a droite et a gauche de la maille 60 61 ijb=ijb_x 62 ije=ije_x 63 64 if (pole_nord.and.ijb==1) ijb=ijb+iip1 65 if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1 66 67 IF (pente_max.gt.-1.e-5) THEN 68 ! IF (pente_max.gt.10) THEN 69 70 ! calcul des pentes avec limitation, Van Leer scheme I: 71 ! ----------------------------------------------------- 72 ! ! on a besoin de q entre ijb et ije 73 ! calcul de la pente aux points u 74 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 75 DO l = 1, llm 76 77 DO ij=ijb,ije-1 78 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 79 ! IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0' 80 ! sigu(ij)=u_m(ij,l)/masse(ij,l,iq) 81 ENDDO 82 DO ij=ijb+iip1-1,ije,iip1 83 dxqu(ij)=dxqu(ij-iim) 84 ! sigu(ij)=sigu(ij-iim) 85 ENDDO 86 87 DO ij=ijb,ije 88 adxqu(ij)=abs(dxqu(ij)) 89 ENDDO 90 91 ! calcul de la pente maximum dans la maille en valeur absolue 92 93 DO ij=ijb+1,ije 94 dxqmax(ij,l)=pente_max* & 95 min(adxqu(ij-1),adxqu(ij)) 96 ! limitation subtile 97 ! , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij))) 98 99 100 ENDDO 101 102 DO ij=ijb+iip1-1,ije,iip1 103 dxqmax(ij-iim,l)=dxqmax(ij,l) 104 ENDDO 105 106 DO ij=ijb+1,ije 107 107 #ifdef CRAY 108 dxq(ij,l)=109 ,cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))108 dxq(ij,l)= & 109 cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij)) 110 110 #else 111 112 113 114 cextremum local115 116 111 IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN 112 dxq(ij,l)=dxqu(ij-1)+dxqu(ij) 113 ELSE 114 ! extremum local 115 dxq(ij,l)=0. 116 ENDIF 117 117 #endif 118 119 dxq(ij,l)=120 ,sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))121 122 123 124 c$OMP END DO NOWAIT125 cprint*,'Ok calcul des pentes'126 127 128 129 cPentes produits:130 c----------------131 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 cextremum local147 148 149 150 151 152 c$OMP END DO NOWAIT153 154 155 156 157 cbouclage de la pente en iip1:158 c-----------------------------159 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)160 161 162 163 164 165 166 167 168 169 c$OMP END DO NOWAIT170 cprint*,'Bouclage en iip1'171 172 ccalcul des flux a gauche et a droite118 dxq(ij,l)=0.5*dxq(ij,l) 119 dxq(ij,l)= & 120 sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l)) 121 ENDDO 122 123 ENDDO ! l=1,llm 124 !$OMP END DO NOWAIT 125 ! print*,'Ok calcul des pentes' 126 127 ELSE ! (pente_max.lt.-1.e-5) 128 129 ! Pentes produits: 130 ! ---------------- 131 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 132 DO l = 1, llm 133 DO ij=ijb,ije-1 134 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 135 ENDDO 136 DO ij=ijb+iip1-1,ije,iip1 137 dxqu(ij)=dxqu(ij-iim) 138 ENDDO 139 140 DO ij=ijb+1,ije 141 zz(ij)=dxqu(ij-1)*dxqu(ij) 142 zz(ij)=zz(ij)+zz(ij) 143 IF(zz(ij).gt.0) THEN 144 dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij)) 145 ELSE 146 ! extremum local 147 dxq(ij,l)=0. 148 ENDIF 149 ENDDO 150 151 ENDDO 152 !$OMP END DO NOWAIT 153 ENDIF ! (pente_max.lt.-1.e-5) 154 155 ! !write(*,*) 'vlx 156: iq,ijb_x=',iq,ijb_x 156 157 ! bouclage de la pente en iip1: 158 ! ----------------------------- 159 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 160 DO l=1,llm 161 DO ij=ijb+iip1-1,ije,iip1 162 dxq(ij-iim,l)=dxq(ij,l) 163 ENDDO 164 DO ij=ijb,ije 165 iadvplus(ij,l)=0 166 ENDDO 167 168 ENDDO 169 !$OMP END DO NOWAIT 170 ! print*,'Bouclage en iip1' 171 172 ! calcul des flux a gauche et a droite 173 173 174 174 #ifdef CRAY 175 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)176 177 178 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),179 , 1.+u_m(ij,l)/masse(ij+1,l,iq),180 ,u_m(ij,l,iq))181 182 u_mq(ij,l)=cvmgp(183 , q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),184 , q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),185 ,u_m(ij,l))186 187 188 189 c$OMP END DO NOWAIT175 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 176 DO l=1,llm 177 DO ij=ijb,ije-1 178 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq), & 179 1.+u_m(ij,l)/masse(ij+1,l,iq), & 180 u_m(ij,l,iq)) 181 zdum(ij,l)=0.5*zdum(ij,l) 182 u_mq(ij,l)=cvmgp( & 183 q(ij,l,iq)+zdum(ij,l)*dxq(ij,l), & 184 q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l), & 185 u_m(ij,l)) 186 u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l) 187 ENDDO 188 ENDDO 189 !$OMP END DO NOWAIT 190 190 #else 191 con cumule le flux correspondant a toutes les mailles dont la masse192 cau travers de la paroi pENDant le pas de temps.193 cprint*,'Cumule ....'194 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)195 196 197 198 cprint*,'masse(',ij,')=',masse(ij,l,iq)199 200 201 u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)202 :+0.5*zdum(ij,l)*dxq(ij,l))203 204 205 u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq)206 :-0.5*zdum(ij,l)*dxq(ij+1,l))207 208 209 210 c$OMP END DO NOWAIT191 ! on cumule le flux correspondant a toutes les mailles dont la masse 192 ! au travers de la paroi pENDant le pas de temps. 193 ! print*,'Cumule ....' 194 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 195 ! ! on a besoin de masse entre ijb et ije 196 DO l=1,llm 197 DO ij=ijb,ije-1 198 ! print*,'masse(',ij,')=',masse(ij,l,iq) 199 IF (u_m(ij,l).gt.0.) THEN 200 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq) 201 u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq) & 202 +0.5*zdum(ij,l)*dxq(ij,l)) 203 ELSE 204 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq) 205 u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq) & 206 -0.5*zdum(ij,l)*dxq(ij+1,l)) 207 ENDIF 208 ENDDO 209 ENDDO 210 !$OMP END DO NOWAIT 211 211 #endif 212 212 213 c go to 9999 214 c detection des points ou on advecte plus que la masse de la 215 c maille 216 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 217 DO l=1,llm 218 DO ij=ijb,ije-1 219 IF(zdum(ij,l).lt.0) THEN 220 iadvplus(ij,l)=1 221 u_mq(ij,l)=0. 222 ENDIF 223 ENDDO 224 ENDDO 225 c$OMP END DO NOWAIT 226 c print*,'Ok test 1' 227 228 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 229 DO l=1,llm 230 DO ij=ijb+iip1-1,ije,iip1 231 iadvplus(ij,l)=iadvplus(ij-iim,l) 232 ENDDO 233 ENDDO 234 c$OMP END DO NOWAIT 235 c print*,'Ok test 2' 236 237 238 c traitement special pour le cas ou on advecte en longitude plus que le 239 c contenu de la maille. 240 c cette partie est mal vectorisee. 241 242 c calcul du nombre de maille sur lequel on advecte plus que la maille. 243 244 n0=0 245 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 246 DO l=1,llm 247 nl(l)=0 248 DO ij=ijb,ije 249 nl(l)=nl(l)+iadvplus(ij,l) 250 ENDDO 251 n0=n0+nl(l) 252 ENDDO 253 c$OMP END DO NOWAIT 254 cym IF(n0.gt.1) THEN 255 cym IF(n0.gt.0) THEN 256 257 c PRINT*,'Nombre de points pour lesquels on advect plus que le' 258 c & ,'contenu de la maille : ',n0 259 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 260 261 262 DO l=1,llm 263 IF(nl(l).gt.0) THEN 264 iju=0 265 c indicage des mailles concernees par le traitement special 266 DO ij=ijb,ije 267 IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN 268 iju=iju+1 269 indu(iju)=ij 270 ENDIF 271 ENDDO 272 niju=iju 273 !PRINT*,'vlx 278, niju,nl',niju,nl(l) 274 275 c traitement des mailles 276 DO iju=1,niju 277 ij=indu(iju) 278 j=(ij-1)/iip1+1 279 zu_m=u_m(ij,l) 280 u_mq(ij,l)=0. 281 IF(zu_m.gt.0.) THEN 282 ijq=ij 283 i=ijq-(j-1)*iip1 284 c accumulation pour les mailles completements advectees 285 do while(zu_m.gt.masse(ijq,l,iq)) 286 u_mq(ij,l)=u_mq(ij,l) 287 & +q(ijq,l,iq)*masse(ijq,l,iq) 288 zu_m=zu_m-masse(ijq,l,iq) 289 i=mod(i-2+iim,iim)+1 290 ijq=(j-1)*iip1+i 291 ENDDO 292 c ajout de la maille non completement advectee 293 u_mq(ij,l)=u_mq(ij,l)+zu_m* 294 & (q(ijq,l,iq)+0.5* 295 & (1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 296 ELSE 297 ijq=ij+1 298 i=ijq-(j-1)*iip1 299 c accumulation pour les mailles completements advectees 300 do while(-zu_m.gt.masse(ijq,l,iq)) 301 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq) 302 & *masse(ijq,l,iq) 303 zu_m=zu_m+masse(ijq,l,iq) 304 i=mod(i,iim)+1 305 ijq=(j-1)*iip1+i 306 ENDDO 307 c ajout de la maille non completement advectee 308 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)- 309 & 0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 310 ENDIF 311 ENDDO 312 ENDIF 313 ENDDO 314 c$OMP END DO NOWAIT 315 cym ENDIF ! n0.gt.0 316 9999 continue 317 318 c bouclage en latitude 319 c print*,'Avant bouclage en latitude' 320 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 321 DO l=1,llm 322 DO ij=ijb+iip1-1,ije,iip1 323 u_mq(ij,l)=u_mq(ij-iim,l) 324 ENDDO 325 ENDDO 326 c$OMP END DO NOWAIT 327 328 ! CRisi: appel récursif de l'advection sur les fils. 329 ! Il faut faire ça avant d'avoir mis à jour q et masse 330 331 do ifils=1,tracers(iq)%nqDescen 332 ! attention: comme Ratio est utilisé comme q dans l'appel 333 ! recursif, il doit contenir à lui seul tous les indices de tous 334 ! les descendants! 335 iq2=tracers(iq)%iqDescen(ifils) 336 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 337 DO l=1,llm 338 DO ij=ijb,ije 339 ! On a besoin de q et masse seulement entre ijb et ije. On ne 340 ! les calcule donc que de ijb à ije 341 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 342 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass) 343 if (q(ij,l,iq).gt.min_qParent) then ! modif 13 nov 2020 344 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 345 else 346 Ratio(ij,l,iq2)=min_ratio 347 endif 348 enddo 349 enddo 350 c$OMP END DO NOWAIT 351 enddo !do ifils=1,tracers(iq)%nqDescen 352 do ifils=1,tracers(iq)%nqChildren 353 iq2=tracers(iq)%iqDescen(ifils) 354 call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2) 213 ! go to 9999 214 ! detection des points ou on advecte plus que la masse de la 215 ! maille 216 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 217 DO l=1,llm 218 DO ij=ijb,ije-1 219 IF(zdum(ij,l).lt.0) THEN 220 iadvplus(ij,l)=1 221 u_mq(ij,l)=0. 222 ENDIF 223 ENDDO 224 ENDDO 225 !$OMP END DO NOWAIT 226 ! print*,'Ok test 1' 227 228 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 229 DO l=1,llm 230 DO ij=ijb+iip1-1,ije,iip1 231 iadvplus(ij,l)=iadvplus(ij-iim,l) 232 ENDDO 233 ENDDO 234 !$OMP END DO NOWAIT 235 ! print*,'Ok test 2' 236 237 238 ! traitement special pour le cas ou on advecte en longitude plus que le 239 ! contenu de la maille. 240 ! cette partie est mal vectorisee. 241 242 ! calcul du nombre de maille sur lequel on advecte plus que la maille. 243 244 n0=0 245 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 246 DO l=1,llm 247 nl(l)=0 248 DO ij=ijb,ije 249 nl(l)=nl(l)+iadvplus(ij,l) 250 ENDDO 251 n0=n0+nl(l) 252 ENDDO 253 !$OMP END DO NOWAIT 254 !ym IF(n0.gt.1) THEN 255 !ym IF(n0.gt.0) THEN 256 257 ! PRINT*,'Nombre de points pour lesquels on advect plus que le' 258 ! & ,'contenu de la maille : ',n0 259 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 260 261 262 DO l=1,llm 263 IF(nl(l).gt.0) THEN 264 iju=0 265 ! indicage des mailles concernees par le traitement special 266 DO ij=ijb,ije 267 IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN 268 iju=iju+1 269 indu(iju)=ij 270 ENDIF 271 ENDDO 272 niju=iju 273 ! !PRINT*,'vlx 278, niju,nl',niju,nl(l) 274 275 ! traitement des mailles 276 DO iju=1,niju 277 ij=indu(iju) 278 j=(ij-1)/iip1+1 279 zu_m=u_m(ij,l) 280 u_mq(ij,l)=0. 281 IF(zu_m.gt.0.) THEN 282 ijq=ij 283 i=ijq-(j-1)*iip1 284 ! accumulation pour les mailles completements advectees 285 do while(zu_m.gt.masse(ijq,l,iq)) 286 u_mq(ij,l)=u_mq(ij,l) & 287 +q(ijq,l,iq)*masse(ijq,l,iq) 288 zu_m=zu_m-masse(ijq,l,iq) 289 i=mod(i-2+iim,iim)+1 290 ijq=(j-1)*iip1+i 291 ENDDO 292 ! ajout de la maille non completement advectee 293 u_mq(ij,l)=u_mq(ij,l)+zu_m* & 294 (q(ijq,l,iq)+0.5* & 295 (1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 296 ELSE 297 ijq=ij+1 298 i=ijq-(j-1)*iip1 299 ! accumulation pour les mailles completements advectees 300 do while(-zu_m.gt.masse(ijq,l,iq)) 301 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq) & 302 *masse(ijq,l,iq) 303 zu_m=zu_m+masse(ijq,l,iq) 304 i=mod(i,iim)+1 305 ijq=(j-1)*iip1+i 306 ENDDO 307 ! ajout de la maille non completement advectee 308 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)- & 309 0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 310 ENDIF 311 ENDDO 312 ENDIF 313 ENDDO 314 !$OMP END DO NOWAIT 315 !ym ENDIF ! n0.gt.0 316 9999 continue 317 318 ! bouclage en latitude 319 ! print*,'Avant bouclage en latitude' 320 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 321 DO l=1,llm 322 DO ij=ijb+iip1-1,ije,iip1 323 u_mq(ij,l)=u_mq(ij-iim,l) 324 ENDDO 325 ENDDO 326 !$OMP END DO NOWAIT 327 328 ! CRisi: appel récursif de l'advection sur les fils. 329 ! Il faut faire ça avant d'avoir mis à jour q et masse 330 331 do ifils=1,tracers(iq)%nqDescen 332 ! ! attention: comme Ratio est utilisé comme q dans l'appel 333 ! ! recursif, il doit contenir à lui seul tous les indices de tous 334 ! ! les descendants! 335 iq2=tracers(iq)%iqDescen(ifils) 336 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 337 DO l=1,llm 338 DO ij=ijb,ije 339 ! ! On a besoin de q et masse seulement entre ijb et ije. On ne 340 ! ! les calcule donc que de ijb à ije 341 ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul 342 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass) 343 if (q(ij,l,iq).gt.min_qParent) then ! modif 13 nov 2020 344 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 345 else 346 Ratio(ij,l,iq2)=min_ratio 347 endif 355 348 enddo 356 ! end CRisi 357 358 359 c calcul des tENDances 360 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 361 DO l=1,llm 362 DO ij=ijb+1,ije 363 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 364 new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),min_qMass) 365 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+ 366 & u_mq(ij-1,l)-u_mq(ij,l)) 367 & /new_m 368 masse(ij,l,iq)=new_m 369 ENDDO 370 c ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous) 371 DO ij=ijb+iip1-1,ije,iip1 372 q(ij-iim,l,iq)=q(ij,l,iq) 373 masse(ij-iim,l,iq)=masse(ij,l,iq) 374 ENDDO 375 ENDDO 376 c$OMP END DO NOWAIT 377 378 ! retablir les fils en rapport de melange par rapport a l'air: 379 ! On calcule q entre ijb+1 et ije -> on fait pareil pour ratio 380 ! puis on boucle en longitude 381 do ifils=1,tracers(iq)%nqDescen 382 iq2=tracers(iq)%iqDescen(ifils) 383 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 384 DO l=1,llm 385 DO ij=ijb+1,ije 386 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 387 enddo 388 DO ij=ijb+iip1-1,ije,iip1 389 q(ij-iim,l,iq2)=q(ij,l,iq2) 390 enddo 391 enddo 392 c$OMP END DO NOWAIT 349 enddo 350 !$OMP END DO NOWAIT 351 enddo !do ifils=1,tracers(iq)%nqDescen 352 do ifils=1,tracers(iq)%nqChildren 353 iq2=tracers(iq)%iqDescen(ifils) 354 call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2) 355 enddo 356 ! end CRisi 357 358 359 ! calcul des tENDances 360 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 361 DO l=1,llm 362 DO ij=ijb+1,ije 363 ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul 364 new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),min_qMass) 365 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+ & 366 u_mq(ij-1,l)-u_mq(ij,l)) & 367 /new_m 368 masse(ij,l,iq)=new_m 369 ENDDO 370 ! ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous) 371 DO ij=ijb+iip1-1,ije,iip1 372 q(ij-iim,l,iq)=q(ij,l,iq) 373 masse(ij-iim,l,iq)=masse(ij,l,iq) 374 ENDDO 375 ENDDO 376 !$OMP END DO NOWAIT 377 378 ! retablir les fils en rapport de melange par rapport a l'air: 379 ! ! On calcule q entre ijb+1 et ije -> on fait pareil pour ratio 380 ! ! puis on boucle en longitude 381 do ifils=1,tracers(iq)%nqDescen 382 iq2=tracers(iq)%iqDescen(ifils) 383 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 384 DO l=1,llm 385 DO ij=ijb+1,ije 386 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 393 387 enddo 394 395 !write(*,*) 'vlsplt 399: iq,ijb_x=',iq,ijb_x 396 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) 397 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1) 398 399 400 RETURN 401 END 402 403 404 RECURSIVE SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v,iq) 405 c 406 c Auteurs: P.Le Van, F.Hourdin, F.Forget 407 c 408 c ******************************************************************** 409 c Shema d'advection " pseudo amont " . 410 c ******************************************************************** 411 c q,masse_adv_v,w sont des arguments d'entree pour le s-pg .... 412 c dq sont des arguments de sortie pour le s-pg .... 413 c 414 c 415 c -------------------------------------------------------------------- 416 USE parallel_lmdz 417 USE infotrac, ONLY : nqtot,tracers, ! CRisi & 418 & min_qParent,min_qMass,min_ratio ! MVals et CRisi 419 USE comconst_mod, ONLY: pi 420 IMPLICIT NONE 421 c 422 include "dimensions.h" 423 include "paramet.h" 424 include "comgeom.h" 425 c 426 c 427 c Arguments: 428 c ---------- 429 REAL masse(ijb_u:ije_u,llm,nqtot),pente_max 430 REAL masse_adv_v( ijb_v:ije_v,llm) 431 REAL q(ijb_u:ije_u,llm,nqtot), dq( ijb_u:ije_u,llm) 432 INTEGER iq ! CRisi 433 c 434 c Local 435 c --------- 436 c 437 INTEGER i,ij,l 438 c 439 REAL airej2,airejjm,airescb(iim),airesch(iim) 440 REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v),zdvm(ijb_u:ije_u,llm) 441 REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u) 442 REAL qbyv(ijb_v:ije_v,llm) 443 444 REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs 445 c REAL newq,oldmasse 446 Logical extremum,first,testcpu 447 REAL temps0,temps1,temps2,temps3,temps4,temps5,second 448 SAVE temps0,temps1,temps2,temps3,temps4,temps5 449 c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5) 450 SAVE first,testcpu 451 c$OMP THREADPRIVATE(first,testcpu) 452 453 REAL convpn,convps,convmpn,convmps 454 real massepn,masseps,qpn,qps 455 REAL sinlon(iip1),sinlondlon(iip1) 456 REAL coslon(iip1),coslondlon(iip1) 457 SAVE sinlon,coslon,sinlondlon,coslondlon 458 c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon) 459 SAVE airej2,airejjm 460 c$OMP THREADPRIVATE(airej2,airejjm) 461 462 REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi 463 INTEGER ifils,iq2 ! CRisi 464 c 465 c 466 REAL SSUM 467 EXTERNAL SSUM 468 469 DATA first,testcpu/.true.,.false./ 470 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./ 471 INTEGER ijb,ije 472 INTEGER ijbm,ijem 473 474 ijb=ij_begin-2*iip1 475 ije=ij_end+2*iip1 476 if (pole_nord) ijb=ij_begin 477 if (pole_sud) ije=ij_end 478 479 IF(first) THEN 480 PRINT*,'Shema Amont nouveau appele dans Vanleer ' 481 first=.false. 482 do i=2,iip1 483 coslon(i)=cos(rlonv(i)) 484 sinlon(i)=sin(rlonv(i)) 485 coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi 486 sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi 487 ENDDO 488 coslon(1)=coslon(iip1) 489 coslondlon(1)=coslondlon(iip1) 490 sinlon(1)=sinlon(iip1) 491 sinlondlon(1)=sinlondlon(iip1) 492 airej2 = SSUM( iim, aire(iip2), 1 ) 493 airejjm= SSUM( iim, aire(ip1jm -iim), 1 ) 388 DO ij=ijb+iip1-1,ije,iip1 389 q(ij-iim,l,iq2)=q(ij,l,iq2) 390 enddo 391 enddo 392 !$OMP END DO NOWAIT 393 enddo 394 395 ! !write(*,*) 'vlsplt 399: iq,ijb_x=',iq,ijb_x 396 ! CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) 397 ! CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1) 398 399 400 RETURN 401 END SUBROUTINE vlx_loc 402 403 404 RECURSIVE SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v,iq) 405 ! 406 ! Auteurs: P.Le Van, F.Hourdin, F.Forget 407 ! 408 ! ******************************************************************** 409 ! Shema d'advection " pseudo amont " . 410 ! ******************************************************************** 411 ! q,masse_adv_v,w sont des arguments d'entree pour le s-pg .... 412 ! dq sont des arguments de sortie pour le s-pg .... 413 ! 414 ! 415 ! -------------------------------------------------------------------- 416 USE parallel_lmdz 417 USE infotrac, ONLY : nqtot,tracers, & ! CRisi & 418 min_qParent,min_qMass,min_ratio ! MVals et CRisi 419 USE comconst_mod, ONLY: pi 420 IMPLICIT NONE 421 ! 422 include "dimensions.h" 423 include "paramet.h" 424 include "comgeom.h" 425 ! 426 ! 427 ! Arguments: 428 ! ---------- 429 REAL :: masse(ijb_u:ije_u,llm,nqtot),pente_max 430 REAL :: masse_adv_v( ijb_v:ije_v,llm) 431 REAL :: q(ijb_u:ije_u,llm,nqtot), dq( ijb_u:ije_u,llm) 432 INTEGER :: iq ! CRisi 433 ! 434 ! Local 435 ! --------- 436 ! 437 INTEGER :: i,ij,l 438 ! 439 REAL :: airej2,airejjm,airescb(iim),airesch(iim) 440 REAL :: dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v),zdvm(ijb_u:ije_u,llm) 441 REAL :: adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u) 442 REAL :: qbyv(ijb_v:ije_v,llm) 443 444 REAL :: qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs 445 ! REAL newq,oldmasse 446 Logical :: extremum,first,testcpu 447 REAL :: temps0,temps1,temps2,temps3,temps4,temps5,second 448 SAVE temps0,temps1,temps2,temps3,temps4,temps5 449 !$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5) 450 SAVE first,testcpu 451 !$OMP THREADPRIVATE(first,testcpu) 452 453 REAL :: convpn,convps,convmpn,convmps 454 real :: massepn,masseps,qpn,qps 455 REAL :: sinlon(iip1),sinlondlon(iip1) 456 REAL :: coslon(iip1),coslondlon(iip1) 457 SAVE sinlon,coslon,sinlondlon,coslondlon 458 !$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon) 459 SAVE airej2,airejjm 460 !$OMP THREADPRIVATE(airej2,airejjm) 461 462 REAL :: Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi 463 INTEGER :: ifils,iq2 ! CRisi 464 ! 465 ! 466 REAL :: SSUM 467 EXTERNAL SSUM 468 469 DATA first,testcpu/.true.,.false./ 470 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./ 471 INTEGER :: ijb,ije 472 INTEGER :: ijbm,ijem 473 474 ijb=ij_begin-2*iip1 475 ije=ij_end+2*iip1 476 if (pole_nord) ijb=ij_begin 477 if (pole_sud) ije=ij_end 478 479 IF(first) THEN 480 PRINT*,'Shema Amont nouveau appele dans Vanleer ' 481 first=.false. 482 do i=2,iip1 483 coslon(i)=cos(rlonv(i)) 484 sinlon(i)=sin(rlonv(i)) 485 coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi 486 sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi 487 ENDDO 488 coslon(1)=coslon(iip1) 489 coslondlon(1)=coslondlon(iip1) 490 sinlon(1)=sinlon(iip1) 491 sinlondlon(1)=sinlondlon(iip1) 492 airej2 = SSUM( iim, aire(iip2), 1 ) 493 airejjm= SSUM( iim, aire(ip1jm -iim), 1 ) 494 ENDIF 495 496 ! 497 ! PRINT*,'CALCUL EN LATITUDE' 498 499 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 500 DO l = 1, llm 501 ! 502 ! -------------------------------- 503 ! CALCUL EN LATITUDE 504 ! -------------------------------- 505 506 ! On commence par calculer la valeur du traceur moyenne sur le premier cercle 507 ! de latitude autour du pole (qpns pour le pole nord et qpsn pour 508 ! le pole nord) qui sera utilisee pour evaluer les pentes au pole. 509 510 if (pole_nord) then 511 DO i = 1, iim 512 airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq) 513 ENDDO 514 qpns = SSUM( iim, airescb ,1 ) / airej2 515 endif 516 517 if (pole_sud) then 518 DO i = 1, iim 519 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq) 520 ENDDO 521 qpsn = SSUM( iim, airesch ,1 ) / airejjm 522 endif 523 524 ! calcul des pentes aux points v 525 526 ijb=ij_begin-2*iip1 527 ije=ij_end+iip1 528 if (pole_nord) ijb=ij_begin 529 if (pole_sud) ije=ij_end-iip1 530 531 ! ! on a besoin de q entre ij_begin-2*iip1 et ij_end+2*iip1 532 ! ! Si pole sud, entre ij_begin-2*iip1 et ij_end 533 ! ! Si pole Nord, entre ij_begin et ij_end+2*iip1 534 DO ij=ijb,ije 535 dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq) 536 adyqv(ij)=abs(dyqv(ij)) 537 ENDDO 538 539 540 ! calcul des pentes aux points scalaires 541 ijb=ij_begin-iip1 542 ije=ij_end+iip1 543 if (pole_nord) ijb=ij_begin+iip1 544 if (pole_sud) ije=ij_end-iip1 545 546 DO ij=ijb,ije 547 dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij)) 548 dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij)) 549 dyqmax(ij)=pente_max*dyqmax(ij) 550 ENDDO 551 552 ! calcul des pentes aux poles 553 IF (pole_nord) THEN 554 DO ij=1,iip1 555 dyq(ij,l)=qpns-q(ij+iip1,l,iq) 556 ENDDO 557 558 dyn1=0. 559 dyn2=0. 560 DO ij=1,iim 561 dyn1=dyn1+sinlondlon(ij)*dyq(ij,l) 562 dyn2=dyn2+coslondlon(ij)*dyq(ij,l) 563 ENDDO 564 DO ij=1,iip1 565 dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij) 566 ENDDO 567 568 DO ij=1,iip1 569 dyq(ij,l)=0. 570 ENDDO 571 ! ym tout cela ne sert pas a grand chose 572 ENDIF 573 574 IF (pole_sud) THEN 575 576 DO ij=1,iip1 577 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn 578 ENDDO 579 580 dys1=0. 581 dys2=0. 582 583 DO ij=1,iim 584 dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l) 585 dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l) 586 ENDDO 587 588 DO ij=1,iip1 589 dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij) 590 ENDDO 591 592 DO ij=1,iip1 593 dyq(ip1jm+ij,l)=0. 594 ENDDO 595 ! ym tout cela ne sert pas a grand chose 596 ENDIF 597 598 ! filtrage de la derivee 599 600 ! calcul des pentes limites aux poles 601 ! ym partie inutile 602 ! goto 8888 603 ! fn=1. 604 ! fs=1. 605 ! DO ij=1,iim 606 ! IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN 607 ! fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn) 608 ! ENDIF 609 ! IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN 610 ! fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs) 611 ! ENDIF 612 ! ENDDO 613 ! DO ij=1,iip1 614 ! dyq(ij,l)=fn*dyq(ij,l) 615 ! dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l) 616 ! ENDDO 617 ! 8888 continue 618 619 620 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 621 ! En memoire de dIFferents tests sur la 622 ! limitation des pentes aux poles. 623 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 624 ! PRINT*,dyq(1) 625 ! PRINT*,dyqv(iip1+1) 626 ! appn=abs(dyq(1)/dyqv(iip1+1)) 627 ! PRINT*,dyq(ip1jm+1) 628 ! PRINT*,dyqv(ip1jm-iip1+1) 629 ! apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) 630 ! DO ij=2,iim 631 ! appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) 632 ! apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) 633 ! ENDDO 634 ! appn=min(pente_max/appn,1.) 635 ! apps=min(pente_max/apps,1.) 636 ! 637 ! 638 ! cas ou on a un extremum au pole 639 ! 640 ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) 641 ! & appn=0. 642 ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* 643 ! & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) 644 ! & apps=0. 645 ! 646 ! limitation des pentes aux poles 647 ! DO ij=1,iip1 648 ! dyq(ij)=appn*dyq(ij) 649 ! dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) 650 ! ENDDO 651 ! 652 ! test 653 ! DO ij=1,iip1 654 ! dyq(iip1+ij)=0. 655 ! dyq(ip1jm+ij-iip1)=0. 656 ! ENDDO 657 ! DO ij=1,ip1jmp1 658 ! dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1)) 659 ! ENDDO 660 ! 661 ! changement 10 07 96 662 ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) 663 ! & THEN 664 ! DO ij=1,iip1 665 ! dyqmax(ij)=0. 666 ! ENDDO 667 ! ELSE 668 ! DO ij=1,iip1 669 ! dyqmax(ij)=pente_max*abs(dyqv(ij)) 670 ! ENDDO 671 ! ENDIF 672 ! 673 ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* 674 ! & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) 675 ! &THEN 676 ! DO ij=ip1jm+1,ip1jmp1 677 ! dyqmax(ij)=0. 678 ! ENDDO 679 ! ELSE 680 ! DO ij=ip1jm+1,ip1jmp1 681 ! dyqmax(ij)=pente_max*abs(dyqv(ij-iip1)) 682 ! ENDDO 683 ! ENDIF 684 ! fin changement 10 07 96 685 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 686 687 ! calcul des pentes limitees 688 ijb=ij_begin-iip1 689 ije=ij_end+iip1 690 if (pole_nord) ijb=ij_begin+iip1 691 if (pole_sud) ije=ij_end-iip1 692 693 DO ij=ijb,ije 694 IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN 695 dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l)) 696 ELSE 697 dyq(ij,l)=0. 698 ENDIF 699 ENDDO 700 701 ENDDO 702 !$OMP END DO NOWAIT 703 704 ijb=ij_begin-iip1 705 ije=ij_end 706 if (pole_nord) ijb=ij_begin 707 if (pole_sud) ije=ij_end-iip1 708 709 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 710 DO l=1,llm 711 DO ij=ijb,ije 712 IF(masse_adv_v(ij,l).gt.0) THEN 713 qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)* & 714 0.5*(1.-masse_adv_v(ij,l) & 715 /masse(ij+iip1,l,iq)) 716 ELSE 717 qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)* & 718 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) 494 719 ENDIF 495 496 c 497 c PRINT*,'CALCUL EN LATITUDE' 498 499 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 500 DO l = 1, llm 501 c 502 c -------------------------------- 503 c CALCUL EN LATITUDE 504 c -------------------------------- 505 506 c On commence par calculer la valeur du traceur moyenne sur le premier cercle 507 c de latitude autour du pole (qpns pour le pole nord et qpsn pour 508 c le pole nord) qui sera utilisee pour evaluer les pentes au pole. 509 510 if (pole_nord) then 511 DO i = 1, iim 512 airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq) 513 ENDDO 514 qpns = SSUM( iim, airescb ,1 ) / airej2 515 endif 516 517 if (pole_sud) then 518 DO i = 1, iim 519 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq) 520 ENDDO 521 qpsn = SSUM( iim, airesch ,1 ) / airejjm 522 endif 523 524 c calcul des pentes aux points v 525 526 ijb=ij_begin-2*iip1 527 ije=ij_end+iip1 528 if (pole_nord) ijb=ij_begin 529 if (pole_sud) ije=ij_end-iip1 530 531 ! on a besoin de q entre ij_begin-2*iip1 et ij_end+2*iip1 532 ! Si pole sud, entre ij_begin-2*iip1 et ij_end 533 ! Si pole Nord, entre ij_begin et ij_end+2*iip1 720 qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l) 721 ENDDO 722 ENDDO 723 !$OMP END DO NOWAIT 724 725 ! CRisi: appel récursif de l'advection sur les fils. 726 ! Il faut faire ça avant d'avoir mis à jour q et masse 727 ! write(*,*)'vly 689: iq,nqChildren(iq)=',iq,tracers(iq)%nqChildren 728 729 ijb=ij_begin-2*iip1 730 ije=ij_end+2*iip1 731 ijbm=ij_begin-iip1 732 ijem=ij_end+iip1 733 if (pole_nord) ijb=ij_begin 734 if (pole_sud) ije=ij_end 735 if (pole_nord) ijbm=ij_begin 736 if (pole_sud) ijem=ij_end 737 738 do ifils=1,tracers(iq)%nqDescen 739 iq2=tracers(iq)%iqDescen(ifils) 740 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 741 DO l=1,llm 742 ! ! modif des bornes: CRisi 16 nov 2020 743 ! ! d'abord masse avec bornes corrigées 744 DO ij=ijbm,ijem 745 ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul 746 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass) 747 enddo 748 749 ! ! ensuite Ratio avec anciennes bornes 534 750 DO ij=ijb,ije 535 dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq) 536 adyqv(ij)=abs(dyqv(ij)) 537 ENDDO 538 539 540 c calcul des pentes aux points scalaires 541 ijb=ij_begin-iip1 542 ije=ij_end+iip1 543 if (pole_nord) ijb=ij_begin+iip1 544 if (pole_sud) ije=ij_end-iip1 545 751 ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul 752 if (q(ij,l,iq).gt.min_qParent) then ! modif 13 nov 2020 753 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 754 else 755 Ratio(ij,l,iq2)=min_ratio 756 endif 757 enddo !DO ij=ijbm,ijem 758 enddo !DO l=1,llm 759 !$OMP END DO NOWAIT 760 enddo 761 762 do ifils=1,tracers(iq)%nqChildren 763 iq2=tracers(iq)%iqDescen(ifils) 764 call vly_loc(Ratio,pente_max,masse,qbyv,iq2) 765 enddo 766 ! end CRisi 767 768 ijb=ij_begin 769 ije=ij_end 770 if (pole_nord) ijb=ij_begin+iip1 771 if (pole_sud) ije=ij_end-iip1 772 773 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 774 DO l=1,llm 775 DO ij=ijb,ije 776 newmasse=masse(ij,l,iq) & 777 +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l) 778 779 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l) & 780 -qbyv(ij-iip1,l))/newmasse 781 782 masse(ij,l,iq)=newmasse 783 784 ENDDO 785 786 787 !.-. ancienne version 788 ! convpn=SSUM(iim,qbyv(1,l),1)/apoln 789 ! convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln 790 if (pole_nord) then 791 convpn=SSUM(iim,qbyv(1,l),1) 792 convmpn=ssum(iim,masse_adv_v(1,l),1) 793 massepn=ssum(iim,masse(1,l,iq),1) 794 qpn=0. 795 do ij=1,iim 796 qpn=qpn+masse(ij,l,iq)*q(ij,l,iq) 797 enddo 798 qpn=(qpn+convpn)/(massepn+convmpn) 799 do ij=1,iip1 800 q(ij,l,iq)=qpn 801 enddo 802 endif 803 804 ! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols 805 ! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols 806 807 if (pole_sud) then 808 809 convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) 810 convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1) 811 masseps=ssum(iim, masse(ip1jm+1,l,iq),1) 812 qps=0. 813 do ij = ip1jm+1,ip1jmp1-1 814 qps=qps+masse(ij,l,iq)*q(ij,l,iq) 815 enddo 816 qps=(qps+convps)/(masseps+convmps) 817 do ij=ip1jm+1,ip1jmp1 818 q(ij,l,iq)=qps 819 enddo 820 endif 821 !.-. fin ancienne version 822 823 !._. nouvelle version 824 ! convpn=SSUM(iim,qbyv(1,l),1) 825 ! convmpn=ssum(iim,masse_adv_v(1,l),1) 826 ! oldmasse=ssum(iim,masse(1,l),1) 827 ! newmasse=oldmasse+convmpn 828 ! newq=(q(1,l)*oldmasse+convpn)/newmasse 829 ! newmasse=newmasse/apoln 830 ! DO ij = 1,iip1 831 ! q(ij,l)=newq 832 ! masse(ij,l,iq)=newmasse*aire(ij) 833 ! ENDDO 834 ! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) 835 ! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1) 836 ! oldmasse=ssum(iim,masse(ip1jm-iim,l),1) 837 ! newmasse=oldmasse+convmps 838 ! newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse 839 ! newmasse=newmasse/apols 840 ! DO ij = ip1jm+1,ip1jmp1 841 ! q(ij,l)=newq 842 ! masse(ij,l,iq)=newmasse*aire(ij) 843 ! ENDDO 844 !._. fin nouvelle version 845 ENDDO 846 !$OMP END DO NOWAIT 847 848 ! retablir les fils en rapport de melange par rapport a l'air: 849 ijb=ij_begin 850 ije=ij_end 851 ! if (pole_nord) ijb=ij_begin 852 ! if (pole_sud) ije=ij_end 853 854 do ifils=1,tracers(iq)%nqDescen 855 iq2=tracers(iq)%iqDescen(ifils) 856 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 857 DO l=1,llm 546 858 DO ij=ijb,ije 547 dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij)) 548 dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij)) 549 dyqmax(ij)=pente_max*dyqmax(ij) 550 ENDDO 551 552 c calcul des pentes aux poles 553 IF (pole_nord) THEN 554 DO ij=1,iip1 555 dyq(ij,l)=qpns-q(ij+iip1,l,iq) 556 ENDDO 557 558 dyn1=0. 559 dyn2=0. 560 DO ij=1,iim 561 dyn1=dyn1+sinlondlon(ij)*dyq(ij,l) 562 dyn2=dyn2+coslondlon(ij)*dyq(ij,l) 563 ENDDO 564 DO ij=1,iip1 565 dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij) 566 ENDDO 567 568 DO ij=1,iip1 569 dyq(ij,l)=0. 570 ENDDO 571 c ym tout cela ne sert pas a grand chose 859 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 860 enddo 861 enddo 862 !$OMP END DO NOWAIT 863 enddo 864 865 866 RETURN 867 END SUBROUTINE vly_loc 868 869 870 871 RECURSIVE SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x,iq) 872 ! 873 ! Auteurs: P.Le Van, F.Hourdin, F.Forget 874 ! 875 ! ******************************************************************** 876 ! Shema d'advection " pseudo amont " . 877 ! ******************************************************************** 878 ! q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... 879 ! dq sont des arguments de sortie pour le s-pg .... 880 ! 881 ! 882 ! -------------------------------------------------------------------- 883 USE parallel_lmdz 884 USE vlz_mod 885 USE infotrac, ONLY : nqtot,tracers, & ! CRisi & 886 min_qParent,min_qMass,min_ratio ! MVals et CRisi 887 888 IMPLICIT NONE 889 ! 890 include "dimensions.h" 891 include "paramet.h" 892 include "iniprint.h" 893 ! 894 ! 895 ! Arguments: 896 ! ---------- 897 REAL :: masse(ijb_u:ije_u,llm,nqtot),pente_max 898 REAL :: q(ijb_u:ije_u,llm,nqtot) 899 REAL :: w(ijb_u:ije_u,llm+1,nqtot) 900 INTEGER :: iq 901 ! 902 ! Local 903 ! --------- 904 ! 905 INTEGER :: i,ij,l,j,ii 906 907 REAL,DIMENSION(ijb_u:ije_u,llm+1) :: wresi,morig,qorig,dzqorig 908 INTEGER,DIMENSION(ijb_u:ije_u,llm+1) :: lorig 909 INTEGER,SAVE :: countcfl 910 !$OMP THREADPRIVATE(countcfl) 911 ! 912 REAL :: newmasse 913 914 REAL :: dzqmax 915 REAL :: sigw 916 917 LOGICAL :: testcpu 918 SAVE testcpu 919 !$OMP THREADPRIVATE(testcpu) 920 REAL :: temps0,temps1,temps2,temps3,temps4,temps5,second 921 SAVE temps0,temps1,temps2,temps3,temps4,temps5 922 !$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5) 923 924 REAL :: SSUM 925 EXTERNAL SSUM 926 927 DATA testcpu/.false./ 928 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./ 929 INTEGER :: ijb,ije,ijb_x,ije_x 930 LOGICAL,SAVE :: first=.TRUE. 931 !$OMP THREADPRIVATE(first) 932 933 ! !REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi 934 ! ! Ces varibles doivent être déclarées en pointer et en save dans 935 ! ! vlz_loc si on veut qu'elles soient vues par tous les threads. 936 INTEGER :: ifils,iq2 ! CRisi 937 938 939 IF (first) THEN 940 first=.FALSE. 941 ENDIF 942 ! On oriente tout dans le sens de la pression c'est a dire dans le 943 ! sens de W 944 945 ! !write(*,*) 'vlsplt 926: entree dans vlz_loc, iq=',iq 946 #ifdef BIDON 947 IF(testcpu) THEN 948 temps0=second(0.) 949 ENDIF 950 #endif 951 952 ijb=ijb_x 953 ije=ije_x 954 955 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 956 DO l=2,llm 957 DO ij=ijb,ije 958 dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq) 959 adzqw(ij,l)=abs(dzqw(ij,l)) 960 ENDDO 961 ENDDO 962 !$OMP END DO 963 964 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 965 DO l=2,llm-1 966 DO ij=ijb,ije 967 #ifdef CRAY 968 dzq(ij,l)=0.5* & 969 cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1)) 970 #else 971 IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 972 dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1)) 973 ELSE 974 dzq(ij,l)=0. 975 ENDIF 976 #endif 977 dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1)) 978 dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l)) 979 ENDDO 980 ENDDO 981 !$OMP END DO NOWAIT 982 983 !$OMP MASTER 984 DO ij=ijb,ije 985 dzq(ij,1)=0. 986 dzq(ij,llm)=0. 987 ENDDO 988 !$OMP END MASTER 989 !$OMP BARRIER 990 #ifdef BIDON 991 IF(testcpu) THEN 992 temps1=temps1+second(0.)-temps0 993 ENDIF 994 #endif 995 996 !-------------------------------------------------------- 997 ! On repere les points qui violent le CFL (|w| > masse) 998 !-------------------------------------------------------- 999 1000 countcfl=0 1001 ! print*,'vlz nouveau' 1002 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1003 DO l = 2,llm 1004 DO ij = ijb,ije 1005 IF( (w(ij,l,iq)>0.AND.w(ij,l,iq)>masse(ij,l,iq)) & 1006 .OR. (w(ij,l,iq)<=0.AND.ABS(w(ij,l,iq))>masse(ij,l-1,iq)) ) & 1007 countcfl=countcfl+1 1008 ENDDO 1009 ENDDO 1010 !$OMP END DO NOWAIT 1011 1012 ! --------------------------------------------------------------- 1013 ! Identification des mailles ou on viole le CFL : w > masse 1014 ! --------------------------------------------------------------- 1015 1016 IF (countcfl==0) THEN 1017 1018 ! --------------------------------------------------------------- 1019 ! .... calcul des termes d'advection verticale ....... 1020 ! Dans le cas où le |w| < masse partout. 1021 ! Version d'origine 1022 ! Pourrait etre enleve si on voit que le code plus general 1023 ! est aussi rapide 1024 ! --------------------------------------------------------------- 1025 1026 ! calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq 1027 1028 ! !write(*,*) 'vlz 982,ijb,ije=',ijb,ije 1029 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1030 DO l = 1,llm-1 1031 do ij = ijb,ije 1032 IF(w(ij,l+1,iq).gt.0.) THEN 1033 sigw=w(ij,l+1,iq)/masse(ij,l+1,iq) 1034 wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l+1,iq) & 1035 +0.5*(1.-sigw)*dzq(ij,l+1)) 1036 ELSE 1037 sigw=w(ij,l+1,iq)/masse(ij,l,iq) 1038 wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l,iq) & 1039 -0.5*(1.+sigw)*dzq(ij,l)) 572 1040 ENDIF 573 574 IF (pole_sud) THEN 575 576 DO ij=1,iip1 577 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn 578 ENDDO 579 580 dys1=0. 581 dys2=0. 582 583 DO ij=1,iim 584 dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l) 585 dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l) 586 ENDDO 587 588 DO ij=1,iip1 589 dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij) 590 ENDDO 591 592 DO ij=1,iip1 593 dyq(ip1jm+ij,l)=0. 594 ENDDO 595 c ym tout cela ne sert pas a grand chose 1041 ENDDO 1042 ENDDO 1043 !$OMP END DO NOWAIT 1044 ! !write(*,*) 'vlz 1001' 1045 1046 ELSE ! countcfl>=1 1047 1048 IF (prt_level>9) THEN 1049 WRITE(lunout,*)'vlz passage dans le non local' 1050 ENDIF 1051 ! --------------------------------------------------------------- 1052 ! Debut du traitement du cas ou on viole le CFL : w > masse 1053 ! --------------------------------------------------------------- 1054 1055 ! Initialisation 1056 1057 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1058 DO l = 2,llm 1059 DO ij = ijb,ije 1060 wresi(ij,l)=w(ij,l,iq) 1061 wq(ij,l,iq)=0. 1062 IF(w(ij,l,iq).gt.0.) THEN 1063 lorig(ij,l)=l 1064 morig(ij,l)=masse(ij,l,iq) 1065 qorig(ij,l)=q(ij,l,iq) 1066 dzqorig(ij,l)=dzq(ij,l) 1067 ELSE 1068 lorig(ij,l)=l-1 1069 morig(ij,l)=masse(ij,l-1,iq) 1070 qorig(ij,l)=q(ij,l-1,iq) 1071 dzqorig(ij,l)=dzq(ij,l-1) 1072 ENDIF 1073 ENDDO 1074 ENDDO 1075 !$OMP END DO NOWAIT 1076 1077 ! Reindicage vertical en accumulant les flux sur 1078 ! les mailles qui viollent le CFL 1079 ! on itère jusqu'à ce que tous les poins satisfassent 1080 ! le critère 1081 DO WHILE (countcfl>=1) 1082 IF (prt_level>9) THEN 1083 WRITE(lunout,*)'On viole le CFL Vertical sur ',countcfl,' pts' 1084 ENDIF 1085 countcfl=0 1086 1087 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1088 DO l = 2,llm 1089 DO ij = ijb,ije 1090 IF (ABS(wresi(ij,l))>morig(ij,l)) THEN 1091 countcfl=countcfl+1 1092 ! rm : les 8 lignes ci dessous pourraient sans doute s'ecrire 1093 ! avec la fonction sign 1094 IF(w(ij,l,iq)>0.) THEN 1095 wresi(ij,l)=wresi(ij,l)-morig(ij,l) 1096 wq(ij,l,iq)=wq(ij,l,iq)+morig(ij,l)*qorig(ij,l) 1097 lorig(ij,l)=lorig(ij,l)+1 1098 ELSE 1099 wresi(ij,l)=wresi(ij,l)+morig(ij,l) 1100 wq(ij,l,iq)=wq(ij,l,iq)-morig(ij,l)*qorig(ij,l) 1101 lorig(ij,l)=lorig(ij,l)-1 1102 ENDIF 1103 ! ! CRisi 24nov2020: ajout d'un message d'erreur clair au lieu d'un plantage 1104 ! ! pour seg fault 1105 if (lorig(ij,l).eq.0) then 1106 call abort_gcm("vlz in vlsplt_loc", & 1107 "unfixable violation of CFL",1) 1108 endif 1109 morig(ij,l)=masse(ij,lorig(ij,l),iq) 1110 qorig(ij,l)=q(ij,lorig(ij,l),iq) 1111 dzqorig(ij,l)=dzq(ij,lorig(ij,l)) 596 1112 ENDIF 597 598 c filtrage de la derivee 599 600 c calcul des pentes limites aux poles 601 c ym partie inutile 602 c goto 8888 603 c fn=1. 604 c fs=1. 605 c DO ij=1,iim 606 c IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN 607 c fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn) 608 c ENDIF 609 c IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN 610 c fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs) 611 c ENDIF 612 c ENDDO 613 c DO ij=1,iip1 614 c dyq(ij,l)=fn*dyq(ij,l) 615 c dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l) 616 c ENDDO 617 c 8888 continue 618 619 620 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 621 C En memoire de dIFferents tests sur la 622 C limitation des pentes aux poles. 623 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 624 C PRINT*,dyq(1) 625 C PRINT*,dyqv(iip1+1) 626 C appn=abs(dyq(1)/dyqv(iip1+1)) 627 C PRINT*,dyq(ip1jm+1) 628 C PRINT*,dyqv(ip1jm-iip1+1) 629 C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) 630 C DO ij=2,iim 631 C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) 632 C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) 633 C ENDDO 634 C appn=min(pente_max/appn,1.) 635 C apps=min(pente_max/apps,1.) 636 C 637 C 638 C cas ou on a un extremum au pole 639 C 640 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) 641 C & appn=0. 642 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* 643 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) 644 C & apps=0. 645 C 646 C limitation des pentes aux poles 647 C DO ij=1,iip1 648 C dyq(ij)=appn*dyq(ij) 649 C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) 650 C ENDDO 651 C 652 C test 653 C DO ij=1,iip1 654 C dyq(iip1+ij)=0. 655 C dyq(ip1jm+ij-iip1)=0. 656 C ENDDO 657 C DO ij=1,ip1jmp1 658 C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1)) 659 C ENDDO 660 C 661 C changement 10 07 96 662 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) 663 C & THEN 664 C DO ij=1,iip1 665 C dyqmax(ij)=0. 666 C ENDDO 667 C ELSE 668 C DO ij=1,iip1 669 C dyqmax(ij)=pente_max*abs(dyqv(ij)) 670 C ENDDO 671 C ENDIF 672 C 673 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* 674 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) 675 C &THEN 676 C DO ij=ip1jm+1,ip1jmp1 677 C dyqmax(ij)=0. 678 C ENDDO 679 C ELSE 680 C DO ij=ip1jm+1,ip1jmp1 681 C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1)) 682 C ENDDO 683 C ENDIF 684 C fin changement 10 07 96 685 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 686 687 c calcul des pentes limitees 688 ijb=ij_begin-iip1 689 ije=ij_end+iip1 690 if (pole_nord) ijb=ij_begin+iip1 691 if (pole_sud) ije=ij_end-iip1 692 1113 ENDDO 1114 ENDDO 1115 !$OMP END DO NOWAIT 1116 1117 ENDDO ! WHILE (countcfl>=1) 1118 1119 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1120 DO l = 2,llm 1121 do ij = ijb,ije 1122 sigw=wresi(ij,l)/morig(ij,l) 1123 IF(w(ij,l,iq).gt.0.) THEN 1124 wq(ij,l,iq)=wq(ij,l,iq)+wresi(ij,l)*(qorig(ij,l) & 1125 +0.5*(1.-sigw)*dzqorig(ij,l)) 1126 ELSE 1127 wq(ij,l,iq)=wq(ij,l,iq)+wresi(ij,l)*(qorig(ij,l) & 1128 -0.5*(1.+sigw)*dzqorig(ij,l)) 1129 ENDIF 1130 ENDDO 1131 ENDDO 1132 !$OMP END DO NOWAIT 1133 1134 1135 ENDIF ! councfl=0 1136 1137 1138 1139 !$OMP MASTER 1140 DO ij=ijb,ije 1141 wq(ij,llm+1,iq)=0. 1142 wq(ij,1,iq)=0. 1143 ENDDO 1144 !$OMP END MASTER 1145 !$OMP BARRIER 1146 1147 ! CRisi: appel récursif de l'advection sur les fils. 1148 ! Il faut faire ça avant d'avoir mis à jour q et masse 1149 ! write(*,*)'vlsplt 942: iq,nqChildren(iq)=',iq,tracers(iq)%nqChildren 1150 do ifils=1,tracers(iq)%nqDescen 1151 iq2=tracers(iq)%iqDescen(ifils) 1152 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1153 DO l=1,llm 693 1154 DO ij=ijb,ije 694 IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN 695 dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l)) 696 ELSE 697 dyq(ij,l)=0. 698 ENDIF 699 ENDDO 700 701 ENDDO 702 c$OMP END DO NOWAIT 703 704 ijb=ij_begin-iip1 705 ije=ij_end 706 if (pole_nord) ijb=ij_begin 707 if (pole_sud) ije=ij_end-iip1 708 709 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 710 DO l=1,llm 711 DO ij=ijb,ije 712 IF(masse_adv_v(ij,l).gt.0) THEN 713 qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)* 714 , 0.5*(1.-masse_adv_v(ij,l) 715 , /masse(ij+iip1,l,iq)) 716 ELSE 717 qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)* 718 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) 719 ENDIF 720 qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l) 721 ENDDO 722 ENDDO 723 c$OMP END DO NOWAIT 724 725 ! CRisi: appel récursif de l'advection sur les fils. 726 ! Il faut faire ça avant d'avoir mis à jour q et masse 727 ! write(*,*)'vly 689: iq,nqChildren(iq)=',iq,tracers(iq)%nqChildren 728 729 ijb=ij_begin-2*iip1 730 ije=ij_end+2*iip1 731 ijbm=ij_begin-iip1 732 ijem=ij_end+iip1 733 if (pole_nord) ijb=ij_begin 734 if (pole_sud) ije=ij_end 735 if (pole_nord) ijbm=ij_begin 736 if (pole_sud) ijem=ij_end 737 738 do ifils=1,tracers(iq)%nqDescen 739 iq2=tracers(iq)%iqDescen(ifils) 740 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 741 DO l=1,llm 742 ! modif des bornes: CRisi 16 nov 2020 743 ! d'abord masse avec bornes corrigées 744 DO ij=ijbm,ijem 745 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 746 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass) 747 enddo 748 749 ! ensuite Ratio avec anciennes bornes 750 DO ij=ijb,ije 751 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 752 if (q(ij,l,iq).gt.min_qParent) then ! modif 13 nov 2020 753 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 754 else 755 Ratio(ij,l,iq2)=min_ratio 756 endif 757 enddo !DO ij=ijbm,ijem 758 enddo !DO l=1,llm 759 c$OMP END DO NOWAIT 1155 ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul 1156 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass) 1157 if (q(ij,l,iq).gt.min_qParent) then 1158 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 1159 else 1160 Ratio(ij,l,iq2)=min_ratio 1161 endif 1162 ! !wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015 1163 w(ij,l,iq2)=wq(ij,l,iq) 760 1164 enddo 761 762 do ifils=1,tracers(iq)%nqChildren 763 iq2=tracers(iq)%iqDescen(ifils) 764 call vly_loc(Ratio,pente_max,masse,qbyv,iq2) 1165 enddo 1166 !$OMP END DO NOWAIT 1167 enddo 1168 !$OMP BARRIER 1169 1170 do ifils=1,tracers(iq)%nqChildren 1171 iq2=tracers(iq)%iqDescen(ifils) 1172 call vlz_loc(Ratio,pente_max,masse,w,ijb_x,ije_x,iq2) 1173 enddo 1174 ! end CRisi 1175 1176 ! CRisi: On rajoute ici une barrière car on veut être sur que tous les 1177 ! wq soient synchronisés 1178 1179 !$OMP BARRIER 1180 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1181 DO l=1,llm 1182 DO ij=ijb,ije 1183 newmasse=masse(ij,l,iq)+w(ij,l+1,iq)-w(ij,l,iq) 1184 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq) & 1185 +wq(ij,l+1,iq)-wq(ij,l,iq)) & 1186 /newmasse 1187 masse(ij,l,iq)=newmasse 1188 ENDDO 1189 ENDDO 1190 !$OMP END DO NOWAIT 1191 1192 1193 ! retablir les fils en rapport de melange par rapport a l'air: 1194 do ifils=1,tracers(iq)%nqDescen 1195 iq2=tracers(iq)%iqDescen(ifils) 1196 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1197 DO l=1,llm 1198 DO ij=ijb,ije 1199 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 765 1200 enddo 766 ! end CRisi 767 768 ijb=ij_begin 769 ije=ij_end 770 if (pole_nord) ijb=ij_begin+iip1 771 if (pole_sud) ije=ij_end-iip1 772 773 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 774 DO l=1,llm 775 DO ij=ijb,ije 776 newmasse=masse(ij,l,iq) 777 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l) 778 779 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l) 780 & -qbyv(ij-iip1,l))/newmasse 781 782 masse(ij,l,iq)=newmasse 783 784 ENDDO 785 786 787 c.-. ancienne version 788 c convpn=SSUM(iim,qbyv(1,l),1)/apoln 789 c convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln 790 if (pole_nord) then 791 convpn=SSUM(iim,qbyv(1,l),1) 792 convmpn=ssum(iim,masse_adv_v(1,l),1) 793 massepn=ssum(iim,masse(1,l,iq),1) 794 qpn=0. 795 do ij=1,iim 796 qpn=qpn+masse(ij,l,iq)*q(ij,l,iq) 797 enddo 798 qpn=(qpn+convpn)/(massepn+convmpn) 799 do ij=1,iip1 800 q(ij,l,iq)=qpn 801 enddo 802 endif 803 804 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols 805 c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols 806 807 if (pole_sud) then 808 809 convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) 810 convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1) 811 masseps=ssum(iim, masse(ip1jm+1,l,iq),1) 812 qps=0. 813 do ij = ip1jm+1,ip1jmp1-1 814 qps=qps+masse(ij,l,iq)*q(ij,l,iq) 815 enddo 816 qps=(qps+convps)/(masseps+convmps) 817 do ij=ip1jm+1,ip1jmp1 818 q(ij,l,iq)=qps 819 enddo 820 endif 821 c.-. fin ancienne version 822 823 c._. nouvelle version 824 c convpn=SSUM(iim,qbyv(1,l),1) 825 c convmpn=ssum(iim,masse_adv_v(1,l),1) 826 c oldmasse=ssum(iim,masse(1,l),1) 827 c newmasse=oldmasse+convmpn 828 c newq=(q(1,l)*oldmasse+convpn)/newmasse 829 c newmasse=newmasse/apoln 830 c DO ij = 1,iip1 831 c q(ij,l)=newq 832 c masse(ij,l,iq)=newmasse*aire(ij) 833 c ENDDO 834 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) 835 c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1) 836 c oldmasse=ssum(iim,masse(ip1jm-iim,l),1) 837 c newmasse=oldmasse+convmps 838 c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse 839 c newmasse=newmasse/apols 840 c DO ij = ip1jm+1,ip1jmp1 841 c q(ij,l)=newq 842 c masse(ij,l,iq)=newmasse*aire(ij) 843 c ENDDO 844 c._. fin nouvelle version 845 ENDDO 846 c$OMP END DO NOWAIT 847 848 ! retablir les fils en rapport de melange par rapport a l'air: 849 ijb=ij_begin 850 ije=ij_end 851 ! if (pole_nord) ijb=ij_begin 852 ! if (pole_sud) ije=ij_end 853 854 do ifils=1,tracers(iq)%nqDescen 855 iq2=tracers(iq)%iqDescen(ifils) 856 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 857 DO l=1,llm 858 DO ij=ijb,ije 859 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 860 enddo 861 enddo 862 c$OMP END DO NOWAIT 863 enddo 864 865 866 RETURN 867 END 868 869 870 871 RECURSIVE SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x,iq) 872 c 873 c Auteurs: P.Le Van, F.Hourdin, F.Forget 874 c 875 c ******************************************************************** 876 c Shema d'advection " pseudo amont " . 877 c ******************************************************************** 878 c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... 879 c dq sont des arguments de sortie pour le s-pg .... 880 c 881 c 882 c -------------------------------------------------------------------- 883 USE parallel_lmdz 884 USE vlz_mod 885 USE infotrac, ONLY : nqtot,tracers, ! CRisi & 886 & min_qParent,min_qMass,min_ratio ! MVals et CRisi 887 888 IMPLICIT NONE 889 c 890 include "dimensions.h" 891 include "paramet.h" 892 include "iniprint.h" 893 c 894 c 895 c Arguments: 896 c ---------- 897 REAL masse(ijb_u:ije_u,llm,nqtot),pente_max 898 REAL q(ijb_u:ije_u,llm,nqtot) 899 REAL w(ijb_u:ije_u,llm+1,nqtot) 900 INTEGER iq 901 c 902 c Local 903 c --------- 904 c 905 INTEGER i,ij,l,j,ii 906 907 REAL,DIMENSION(ijb_u:ije_u,llm+1) :: wresi,morig,qorig,dzqorig 908 INTEGER,DIMENSION(ijb_u:ije_u,llm+1) :: lorig 909 INTEGER,SAVE :: countcfl 910 !$OMP THREADPRIVATE(countcfl) 911 c 912 REAL newmasse 913 914 REAL dzqmax 915 REAL sigw 916 917 LOGICAL testcpu 918 SAVE testcpu 919 c$OMP THREADPRIVATE(testcpu) 920 REAL temps0,temps1,temps2,temps3,temps4,temps5,second 921 SAVE temps0,temps1,temps2,temps3,temps4,temps5 922 c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5) 923 924 REAL SSUM 925 EXTERNAL SSUM 926 927 DATA testcpu/.false./ 928 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./ 929 INTEGER ijb,ije,ijb_x,ije_x 930 LOGICAL,SAVE :: first=.TRUE. 931 !$OMP THREADPRIVATE(first) 932 933 !REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi 934 ! Ces varibles doivent être déclarées en pointer et en save dans 935 ! vlz_loc si on veut qu'elles soient vues par tous les threads. 936 INTEGER ifils,iq2 ! CRisi 937 938 939 IF (first) THEN 940 first=.FALSE. 941 ENDIF 942 c On oriente tout dans le sens de la pression c'est a dire dans le 943 c sens de W 944 945 !write(*,*) 'vlsplt 926: entree dans vlz_loc, iq=',iq 946 #ifdef BIDON 947 IF(testcpu) THEN 948 temps0=second(0.) 949 ENDIF 950 #endif 951 952 ijb=ijb_x 953 ije=ije_x 954 955 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 956 DO l=2,llm 957 DO ij=ijb,ije 958 dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq) 959 adzqw(ij,l)=abs(dzqw(ij,l)) 960 ENDDO 961 ENDDO 962 c$OMP END DO 963 964 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 965 DO l=2,llm-1 966 DO ij=ijb,ije 967 #ifdef CRAY 968 dzq(ij,l)=0.5* 969 , cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1)) 970 #else 971 IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 972 dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1)) 973 ELSE 974 dzq(ij,l)=0. 975 ENDIF 976 #endif 977 dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1)) 978 dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l)) 979 ENDDO 980 ENDDO 981 c$OMP END DO NOWAIT 982 983 c$OMP MASTER 984 DO ij=ijb,ije 985 dzq(ij,1)=0. 986 dzq(ij,llm)=0. 987 ENDDO 988 c$OMP END MASTER 989 c$OMP BARRIER 990 #ifdef BIDON 991 IF(testcpu) THEN 992 temps1=temps1+second(0.)-temps0 993 ENDIF 994 #endif 995 996 !-------------------------------------------------------- 997 ! On repere les points qui violent le CFL (|w| > masse) 998 !-------------------------------------------------------- 999 1000 countcfl=0 1001 ! print*,'vlz nouveau' 1002 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1003 DO l = 2,llm 1004 DO ij = ijb,ije 1005 IF( (w(ij,l,iq)>0.AND.w(ij,l,iq)>masse(ij,l,iq)) 1006 s .OR. (w(ij,l,iq)<=0.AND.ABS(w(ij,l,iq))>masse(ij,l-1,iq)) ) 1007 s countcfl=countcfl+1 1008 ENDDO 1009 ENDDO 1010 c$OMP END DO NOWAIT 1011 1012 c --------------------------------------------------------------- 1013 c Identification des mailles ou on viole le CFL : w > masse 1014 c --------------------------------------------------------------- 1015 1016 IF (countcfl==0) THEN 1017 1018 c --------------------------------------------------------------- 1019 c .... calcul des termes d'advection verticale ....... 1020 c Dans le cas où le |w| < masse partout. 1021 c Version d'origine 1022 c Pourrait etre enleve si on voit que le code plus general 1023 c est aussi rapide 1024 c --------------------------------------------------------------- 1025 1026 c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq 1027 1028 !write(*,*) 'vlz 982,ijb,ije=',ijb,ije 1029 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1030 DO l = 1,llm-1 1031 do ij = ijb,ije 1032 IF(w(ij,l+1,iq).gt.0.) THEN 1033 sigw=w(ij,l+1,iq)/masse(ij,l+1,iq) 1034 wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l+1,iq) 1035 : +0.5*(1.-sigw)*dzq(ij,l+1)) 1036 ELSE 1037 sigw=w(ij,l+1,iq)/masse(ij,l,iq) 1038 wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l,iq) 1039 : -0.5*(1.+sigw)*dzq(ij,l)) 1040 ENDIF 1041 ENDDO 1042 ENDDO 1043 c$OMP END DO NOWAIT 1044 !write(*,*) 'vlz 1001' 1045 1046 ELSE ! countcfl>=1 1047 1048 IF (prt_level>9) THEN 1049 WRITE(lunout,*)'vlz passage dans le non local' 1050 ENDIF 1051 c --------------------------------------------------------------- 1052 c Debut du traitement du cas ou on viole le CFL : w > masse 1053 c --------------------------------------------------------------- 1054 1055 c Initialisation 1056 1057 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1058 DO l = 2,llm 1059 DO ij = ijb,ije 1060 wresi(ij,l)=w(ij,l,iq) 1061 wq(ij,l,iq)=0. 1062 IF(w(ij,l,iq).gt.0.) THEN 1063 lorig(ij,l)=l 1064 morig(ij,l)=masse(ij,l,iq) 1065 qorig(ij,l)=q(ij,l,iq) 1066 dzqorig(ij,l)=dzq(ij,l) 1067 ELSE 1068 lorig(ij,l)=l-1 1069 morig(ij,l)=masse(ij,l-1,iq) 1070 qorig(ij,l)=q(ij,l-1,iq) 1071 dzqorig(ij,l)=dzq(ij,l-1) 1072 ENDIF 1073 ENDDO 1074 ENDDO 1075 c$OMP END DO NOWAIT 1076 1077 c Reindicage vertical en accumulant les flux sur 1078 c les mailles qui viollent le CFL 1079 c on itère jusqu'à ce que tous les poins satisfassent 1080 c le critère 1081 DO WHILE (countcfl>=1) 1082 IF (prt_level>9) THEN 1083 WRITE(lunout,*)'On viole le CFL Vertical sur ',countcfl,' pts' 1084 ENDIF 1085 countcfl=0 1086 1087 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1088 DO l = 2,llm 1089 DO ij = ijb,ije 1090 IF (ABS(wresi(ij,l))>morig(ij,l)) THEN 1091 countcfl=countcfl+1 1092 ! rm : les 8 lignes ci dessous pourraient sans doute s'ecrire 1093 ! avec la fonction sign 1094 IF(w(ij,l,iq)>0.) THEN 1095 wresi(ij,l)=wresi(ij,l)-morig(ij,l) 1096 wq(ij,l,iq)=wq(ij,l,iq)+morig(ij,l)*qorig(ij,l) 1097 lorig(ij,l)=lorig(ij,l)+1 1098 ELSE 1099 wresi(ij,l)=wresi(ij,l)+morig(ij,l) 1100 wq(ij,l,iq)=wq(ij,l,iq)-morig(ij,l)*qorig(ij,l) 1101 lorig(ij,l)=lorig(ij,l)-1 1102 ENDIF 1103 ! CRisi 24nov2020: ajout d'un message d'erreur clair au lieu d'un plantage 1104 ! pour seg fault 1105 if (lorig(ij,l).eq.0) then 1106 call abort_gcm("vlz in vlsplt_loc", 1107 : "unfixable violation of CFL",1) 1108 endif 1109 morig(ij,l)=masse(ij,lorig(ij,l),iq) 1110 qorig(ij,l)=q(ij,lorig(ij,l),iq) 1111 dzqorig(ij,l)=dzq(ij,lorig(ij,l)) 1112 ENDIF 1113 ENDDO 1114 ENDDO 1115 c$OMP END DO NOWAIT 1116 1117 ENDDO ! WHILE (countcfl>=1) 1118 1119 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1120 DO l = 2,llm 1121 do ij = ijb,ije 1122 sigw=wresi(ij,l)/morig(ij,l) 1123 IF(w(ij,l,iq).gt.0.) THEN 1124 wq(ij,l,iq)=wq(ij,l,iq)+wresi(ij,l)*(qorig(ij,l) 1125 : +0.5*(1.-sigw)*dzqorig(ij,l)) 1126 ELSE 1127 wq(ij,l,iq)=wq(ij,l,iq)+wresi(ij,l)*(qorig(ij,l) 1128 : -0.5*(1.+sigw)*dzqorig(ij,l)) 1129 ENDIF 1130 ENDDO 1131 ENDDO 1132 c$OMP END DO NOWAIT 1133 1134 1135 ENDIF ! councfl=0 1136 1137 1138 1139 c$OMP MASTER 1140 DO ij=ijb,ije 1141 wq(ij,llm+1,iq)=0. 1142 wq(ij,1,iq)=0. 1143 ENDDO 1144 c$OMP END MASTER 1145 c$OMP BARRIER 1146 1147 ! CRisi: appel récursif de l'advection sur les fils. 1148 ! Il faut faire ça avant d'avoir mis à jour q et masse 1149 ! write(*,*)'vlsplt 942: iq,nqChildren(iq)=',iq,tracers(iq)%nqChildren 1150 do ifils=1,tracers(iq)%nqDescen 1151 iq2=tracers(iq)%iqDescen(ifils) 1152 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1153 DO l=1,llm 1154 DO ij=ijb,ije 1155 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 1156 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass) 1157 if (q(ij,l,iq).gt.min_qParent) then 1158 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 1159 else 1160 Ratio(ij,l,iq2)=min_ratio 1161 endif 1162 !wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015 1163 w(ij,l,iq2)=wq(ij,l,iq) 1164 enddo 1165 enddo 1166 c$OMP END DO NOWAIT 1167 enddo 1168 c$OMP BARRIER 1169 1170 do ifils=1,tracers(iq)%nqChildren 1171 iq2=tracers(iq)%iqDescen(ifils) 1172 call vlz_loc(Ratio,pente_max,masse,w,ijb_x,ije_x,iq2) 1173 enddo 1174 ! end CRisi 1175 1176 ! CRisi: On rajoute ici une barrière car on veut être sur que tous les 1177 ! wq soient synchronisés 1178 1179 c$OMP BARRIER 1180 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1181 DO l=1,llm 1182 DO ij=ijb,ije 1183 newmasse=masse(ij,l,iq)+w(ij,l+1,iq)-w(ij,l,iq) 1184 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq) 1185 & +wq(ij,l+1,iq)-wq(ij,l,iq)) 1186 & /newmasse 1187 masse(ij,l,iq)=newmasse 1188 ENDDO 1189 ENDDO 1190 c$OMP END DO NOWAIT 1191 1192 1193 ! retablir les fils en rapport de melange par rapport a l'air: 1194 do ifils=1,tracers(iq)%nqDescen 1195 iq2=tracers(iq)%iqDescen(ifils) 1196 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1197 DO l=1,llm 1198 DO ij=ijb,ije 1199 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 1200 enddo 1201 enddo 1202 c$OMP END DO NOWAIT 1203 enddo 1204 1205 RETURN 1206 END 1207 c SUBROUTINE minmaxq(zq,qmin,qmax,comment) 1208 c 1209 c INCLUDE "dimensions.h" 1210 c INCLUDE "paramet.h" 1211 1212 c CHARACTER*(*) comment 1213 c real qmin,qmax 1214 c real zq(ip1jmp1,llm) 1215 1216 c INTEGER jadrs(ip1jmp1), jbad, k, i 1217 1218 1219 c DO k = 1, llm 1220 c jbad = 0 1221 c DO i = 1, ip1jmp1 1222 c IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN 1223 c jbad = jbad + 1 1224 c jadrs(jbad) = i 1225 c ENDIF 1226 c ENDDO 1227 c IF (jbad.GT.0) THEN 1228 c PRINT*, comment 1229 c DO i = 1, jbad 1230 cc PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k) 1231 c ENDDO 1232 c ENDIF 1233 c ENDDO 1234 1235 c return 1236 c end 1237 1238 1239 1240 1201 enddo 1202 !$OMP END DO NOWAIT 1203 enddo 1204 1205 RETURN 1206 END SUBROUTINE vlz_loc 1207 ! SUBROUTINE minmaxq(zq,qmin,qmax,comment) 1208 ! 1209 ! INCLUDE "dimensions.h" 1210 ! INCLUDE "paramet.h" 1211 1212 ! CHARACTER*(*) comment 1213 ! real qmin,qmax 1214 ! real zq(ip1jmp1,llm) 1215 1216 ! INTEGER jadrs(ip1jmp1), jbad, k, i 1217 1218 1219 ! DO k = 1, llm 1220 ! jbad = 0 1221 ! DO i = 1, ip1jmp1 1222 ! IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN 1223 ! jbad = jbad + 1 1224 ! jadrs(jbad) = i 1225 ! ENDIF 1226 ! ENDDO 1227 ! IF (jbad.GT.0) THEN 1228 ! PRINT*, comment 1229 ! DO i = 1, jbad 1230 !c PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k) 1231 ! ENDDO 1232 ! ENDIF 1233 ! ENDDO 1234 1235 ! return 1236 ! end 1237 1238 1239 1240
Note: See TracChangeset
for help on using the changeset viewer.
