Changeset 5079
- Timestamp:
- Jul 19, 2024, 11:28:59 AM (7 weeks ago)
- Location:
- LMDZ6/trunk/libf
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/dyn3d_common/advy.F
r2600 r5079 121 121 enddo 122 122 123 DO 1L=1,NIV123 DO L=1,NIV 124 124 C 125 125 C place limits on appropriate moments before transport … … 128 128 IF(.NOT.LIMIT) GO TO 11 129 129 C 130 DO 10JV=1,NTRA131 DO 10K=1,LAT132 DO 100I=1,LON130 DO JV=1,NTRA 131 DO K=1,LAT 132 DO I=1,LON 133 133 sy(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.), 134 134 + ABS(sy(I,K,L,JV))),sy(I,K,L,JV)) 135 100 CONTINUE 136 10 CONTINUE 135 END DO 136 END DO 137 END DO 137 138 C 138 139 11 CONTINUE … … 141 142 C 142 143 SM0=0. 143 DO 20JV=1,NTRA144 DO JV=1,NTRA 144 145 S00(JV)=0. 145 20 CONTINUE146 C 147 DO 21I=1,LON148 C 149 IF(VGRI(I,0,L) .LE.0.) THEN146 END DO 147 C 148 DO I=1,LON 149 C 150 IF(VGRI(I,0,L)<=0.) THEN 150 151 FM(I,0)=-VGRI(I,0,L)*DTY 151 152 ALF(I,0)=FM(I,0)/SM(I,1,L) … … 158 159 ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0) 159 160 C 160 21 CONTINUE161 C 162 DO 22JV=1,NTRA163 DO 220I=1,LON164 C 165 IF(VGRI(I,0,L) .LE.0.) THEN161 END DO 162 C 163 DO JV=1,NTRA 164 DO I=1,LON 165 C 166 IF(VGRI(I,0,L)<=0.) THEN 166 167 C 167 168 F0(I,0,JV)=ALF(I,0)* … … 176 177 ENDIF 177 178 C 178 220 CONTINUE179 22 CONTINUE180 C 181 DO 23I=1,LON182 IF(VGRI(I,0,L) .GT.0.) THEN179 END DO 180 END DO 181 C 182 DO I=1,LON 183 IF(VGRI(I,0,L)>0.) THEN 183 184 FM(I,0)=VGRI(I,0,L)*DTY 184 185 ALF(I,0)=FM(I,0)/SM0 185 186 ENDIF 186 23 CONTINUE187 C 188 DO 24JV=1,NTRA189 DO 240I=1,LON190 IF(VGRI(I,0,L) .GT.0.) THEN187 END DO 188 C 189 DO JV=1,NTRA 190 DO I=1,LON 191 IF(VGRI(I,0,L)>0.) THEN 191 192 F0(I,0,JV)=ALF(I,0)*S00(JV) 192 193 ENDIF 193 240 CONTINUE194 24 CONTINUE194 END DO 195 END DO 195 196 C 196 197 C puts the temporary moments Fi into appropriate neighboring boxes 197 198 C 198 DO 25I=1,LON199 C 200 IF(VGRI(I,0,L) .GT.0.) THEN199 DO I=1,LON 200 C 201 IF(VGRI(I,0,L)>0.) THEN 201 202 SM(I,1,L)=SM(I,1,L)+FM(I,0) 202 203 ALF(I,0)=FM(I,0)/SM(I,1,L) … … 205 206 ALF1(I,0)=1.-ALF(I,0) 206 207 C 207 25 CONTINUE208 C 209 DO 26JV=1,NTRA210 DO 260I=1,LON211 C 212 IF(VGRI(I,0,L) .GT.0.) THEN208 END DO 209 C 210 DO JV=1,NTRA 211 DO I=1,LON 212 C 213 IF(VGRI(I,0,L)>0.) THEN 213 214 C 214 215 TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV) … … 218 219 ENDIF 219 220 C 220 260 CONTINUE221 26 CONTINUE221 END DO 222 END DO 222 223 C 223 224 C calculate flux and moments between adjacent boxes … … 227 228 C flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0 228 229 C 229 DO 30K=1,LAT-1230 DO K=1,LAT-1 230 231 KP=K+1 231 DO 300I=1,LON232 C 233 IF(VGRI(I,K,L) .LT.0.) THEN232 DO I=1,LON 233 C 234 IF(VGRI(I,K,L)<0.) THEN 234 235 FM(I,K)=-VGRI(I,K,L)*DTY 235 236 ALF(I,K)=FM(I,K)/SM(I,KP,L) … … 245 246 ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K) 246 247 C 247 300 CONTINUE248 30 CONTINUE249 C 250 DO 31JV=1,NTRA251 DO 31K=1,LAT-1248 END DO 249 END DO 250 C 251 DO JV=1,NTRA 252 DO K=1,LAT-1 252 253 KP=K+1 253 DO 310I=1,LON254 C 255 IF(VGRI(I,K,L) .LT.0.) THEN254 DO I=1,LON 255 C 256 IF(VGRI(I,K,L)<0.) THEN 256 257 C 257 258 F0(I,K,JV)=ALF (I,K)* … … 281 282 ENDIF 282 283 C 283 310 CONTINUE 284 31 CONTINUE 284 END DO 285 END DO 286 END DO 285 287 C 286 288 C puts the temporary moments Fi into appropriate neighboring boxes 287 289 C 288 DO 32K=1,LAT-1290 DO K=1,LAT-1 289 291 KP=K+1 290 DO 320I=1,LON291 C 292 IF(VGRI(I,K,L) .LT.0.) THEN292 DO I=1,LON 293 C 294 IF(VGRI(I,K,L)<0.) THEN 293 295 SM(I,K,L)=SM(I,K,L)+FM(I,K) 294 296 ALF(I,K)=FM(I,K)/SM(I,K,L) … … 300 302 ALF1(I,K)=1.-ALF(I,K) 301 303 C 302 320 CONTINUE303 32 CONTINUE304 C 305 DO 33JV=1,NTRA306 DO 33K=1,LAT-1304 END DO 305 END DO 306 C 307 DO JV=1,NTRA 308 DO K=1,LAT-1 307 309 KP=K+1 308 DO 330I=1,LON309 C 310 IF(VGRI(I,K,L) .LT.0.) THEN310 DO I=1,LON 311 C 312 IF(VGRI(I,K,L)<0.) THEN 311 313 C 312 314 TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV) … … 328 330 ENDIF 329 331 C 330 330 CONTINUE 331 33 CONTINUE 332 END DO 333 END DO 334 END DO 332 335 C 333 336 C traitement special pour le pole Sud (idem pole Nord) … … 336 339 C 337 340 SM0=0. 338 DO 40JV=1,NTRA341 DO JV=1,NTRA 339 342 S00(JV)=0. 340 40 CONTINUE341 C 342 DO 41I=1,LON343 C 344 IF(VGRI(I,K,L) .GE.0.) THEN343 END DO 344 C 345 DO I=1,LON 346 C 347 IF(VGRI(I,K,L)>=0.) THEN 345 348 FM(I,K)=VGRI(I,K,L)*DTY 346 349 ALF(I,K)=FM(I,K)/SM(I,K,L) … … 353 356 ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K) 354 357 C 355 41 CONTINUE356 C 357 DO 42JV=1,NTRA358 DO 420I=1,LON359 C 360 IF(VGRI(I,K,L) .GE.0.) THEN358 END DO 359 C 360 DO JV=1,NTRA 361 DO I=1,LON 362 C 363 IF(VGRI(I,K,L)>=0.) THEN 361 364 F0 (I,K,JV)=ALF(I,K)* 362 365 + ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) ) … … 369 372 ENDIF 370 373 C 371 420 CONTINUE372 42 CONTINUE373 C 374 DO 43I=1,LON375 IF(VGRI(I,K,L) .LT.0.) THEN374 END DO 375 END DO 376 C 377 DO I=1,LON 378 IF(VGRI(I,K,L)<0.) THEN 376 379 FM(I,K)=-VGRI(I,K,L)*DTY 377 380 ALF(I,K)=FM(I,K)/SM0 378 381 ENDIF 379 43 CONTINUE380 C 381 DO 44JV=1,NTRA382 DO 440I=1,LON383 IF(VGRI(I,K,L) .LT.0.) THEN382 END DO 383 C 384 DO JV=1,NTRA 385 DO I=1,LON 386 IF(VGRI(I,K,L)<0.) THEN 384 387 F0(I,K,JV)=ALF(I,K)*S00(JV) 385 388 ENDIF 386 440 CONTINUE387 44 CONTINUE389 END DO 390 END DO 388 391 C 389 392 C puts the temporary moments Fi into appropriate neighboring boxes 390 393 C 391 DO 45I=1,LON392 C 393 IF(VGRI(I,K,L) .LT.0.) THEN394 DO I=1,LON 395 C 396 IF(VGRI(I,K,L)<0.) THEN 394 397 SM(I,K,L)=SM(I,K,L)+FM(I,K) 395 398 ALF(I,K)=FM(I,K)/SM(I,K,L) … … 398 401 ALF1(I,K)=1.-ALF(I,K) 399 402 C 400 45 CONTINUE401 C 402 DO 46JV=1,NTRA403 DO 460I=1,LON404 C 405 IF(VGRI(I,K,L) .LT.0.) THEN403 END DO 404 C 405 DO JV=1,NTRA 406 DO I=1,LON 407 C 408 IF(VGRI(I,K,L)<0.) THEN 406 409 C 407 410 TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV) … … 411 414 ENDIF 412 415 C 413 460 CONTINUE414 46 CONTINUE415 C 416 1 CONTINUE416 END DO 417 END DO 418 C 419 END DO 417 420 C 418 421 RETURN -
LMDZ6/trunk/libf/dyn3d_common/advyp.F
r2600 r5079 153 153 C-AA 20/10/94 le signe -1 est necessaire car indexation opposee 154 154 155 DO 500l = 1,llm156 DO 500j = 1,jjm157 DO 500 i = 1,iip1155 DO l = 1,llm 156 DO j = 1,jjm 157 DO i = 1,iip1 158 158 vgri (i,j,llm+1-l)=-1.*pbarv (i,j,l) 159 500 CONTINUE 159 END DO 160 END DO 161 END DO 160 162 161 163 CAA Initialisation de flux fictifs aux bords sup. des boites pol. … … 171 173 C boucle sur les niveaux 172 174 C 173 DO 1L=1,NIV175 DO L=1,NIV 174 176 C 175 177 C place limits on appropriate moments before transport … … 178 180 IF(.NOT.LIMIT) GO TO 11 179 181 C 180 DO 10JV=1,NTRA181 DO 10K=1,LAT182 DO 100I=1,LON183 IF(S0(I,K,L,JV) .GT.0.) THEN182 DO JV=1,NTRA 183 DO K=1,LAT 184 DO I=1,LON 185 IF(S0(I,K,L,JV)>0.) THEN 184 186 SLPMAX=AMAX1(S0(I,K,L,JV),0.) 185 187 S1MAX=1.5*SLPMAX … … 197 199 SYZ(I,K,L,JV)=0. 198 200 ENDIF 199 100 CONTINUE 200 10 CONTINUE 201 END DO 202 END DO 203 END DO 201 204 C 202 205 11 CONTINUE … … 205 208 C 206 209 SM0=0. 207 DO 20JV=1,NTRA210 DO JV=1,NTRA 208 211 S00(JV)=0. 209 20 CONTINUE210 C 211 DO 21I=1,LON212 C 213 IF(VGRI(I,0,L) .LE.0.) THEN212 END DO 213 C 214 DO I=1,LON 215 C 216 IF(VGRI(I,0,L)<=0.) THEN 214 217 FM(I,0)=-VGRI(I,0,L)*DTY 215 218 ALF(I,0)=FM(I,0)/SM(I,1,L) … … 225 228 ALF4(I,0)=ALF1(I,0)*ALF1Q(I,0) 226 229 C 227 21 CONTINUE230 END DO 228 231 c print*,'ADVYP 21' 229 232 C 230 DO 22JV=1,NTRA231 DO 220I=1,LON232 C 233 IF(VGRI(I,0,L) .LE.0.) THEN233 DO JV=1,NTRA 234 DO I=1,LON 235 C 236 IF(VGRI(I,0,L)<=0.) THEN 234 237 C 235 238 F0(I,0,JV)=ALF(I,0)* ( S0(I,1,L,JV)-ALF1(I,0)* … … 253 256 ENDIF 254 257 C 255 220 CONTINUE256 22 CONTINUE257 C 258 DO 23I=1,LON259 IF(VGRI(I,0,L) .GT.0.) THEN258 END DO 259 END DO 260 C 261 DO I=1,LON 262 IF(VGRI(I,0,L)>0.) THEN 260 263 FM(I,0)=VGRI(I,0,L)*DTY 261 264 ALF(I,0)=FM(I,0)/SM0 262 265 ENDIF 263 23 CONTINUE264 C 265 DO 24JV=1,NTRA266 DO 240I=1,LON267 IF(VGRI(I,0,L) .GT.0.) THEN266 END DO 267 C 268 DO JV=1,NTRA 269 DO I=1,LON 270 IF(VGRI(I,0,L)>0.) THEN 268 271 F0(I,0,JV)=ALF(I,0)*S00(JV) 269 272 ENDIF 270 240 CONTINUE271 24 CONTINUE273 END DO 274 END DO 272 275 C 273 276 C puts the temporary moments Fi into appropriate neighboring boxes 274 277 C 275 278 c print*,'av ADVYP 25' 276 DO 25I=1,LON277 C 278 IF(VGRI(I,0,L) .GT.0.) THEN279 DO I=1,LON 280 C 281 IF(VGRI(I,0,L)>0.) THEN 279 282 SM(I,1,L)=SM(I,1,L)+FM(I,0) 280 283 ALF(I,0)=FM(I,0)/SM(I,1,L) … … 287 290 ALF3(I,0)=ALF1(I,0)*ALF(I,0) 288 291 C 289 25 CONTINUE292 END DO 290 293 c print*,'av ADVYP 25' 291 294 C 292 DO 26JV=1,NTRA293 DO 260I=1,LON294 C 295 IF(VGRI(I,0,L) .GT.0.) THEN295 DO JV=1,NTRA 296 DO I=1,LON 297 C 298 IF(VGRI(I,0,L)>0.) THEN 296 299 C 297 300 TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV) … … 305 308 ENDIF 306 309 C 307 260 CONTINUE308 26 CONTINUE310 END DO 311 END DO 309 312 C 310 313 C calculate flux and moments between adjacent boxes … … 315 318 C 316 319 c print*,'av ADVYP 30' 317 DO 30K=1,LAT-1320 DO K=1,LAT-1 318 321 KP=K+1 319 DO 300I=1,LON320 C 321 IF(VGRI(I,K,L) .LT.0.) THEN322 DO I=1,LON 323 C 324 IF(VGRI(I,K,L)<0.) THEN 322 325 FM(I,K)=-VGRI(I,K,L)*DTY 323 326 ALF(I,K)=FM(I,K)/SM(I,KP,L) … … 336 339 ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K) 337 340 C 338 300 CONTINUE339 30 CONTINUE341 END DO 342 END DO 340 343 c print*,'ap ADVYP 30' 341 344 C 342 DO 31JV=1,NTRA343 DO 31K=1,LAT-1345 DO JV=1,NTRA 346 DO K=1,LAT-1 344 347 KP=K+1 345 DO 310I=1,LON346 C 347 IF(VGRI(I,K,L) .LT.0.) THEN348 DO I=1,LON 349 C 350 IF(VGRI(I,K,L)<0.) THEN 348 351 C 349 352 F0 (I,K,JV)=ALF (I,K)* ( S0(I,KP,L,JV)-ALF1(I,K)* … … 403 406 ENDIF 404 407 C 405 310 CONTINUE 406 31 CONTINUE 408 END DO 409 END DO 410 END DO 407 411 c print*,'ap ADVYP 31' 408 412 C 409 413 C puts the temporary moments Fi into appropriate neighboring boxes 410 414 C 411 DO 32K=1,LAT-1415 DO K=1,LAT-1 412 416 KP=K+1 413 DO 320I=1,LON414 C 415 IF(VGRI(I,K,L) .LT.0.) THEN417 DO I=1,LON 418 C 419 IF(VGRI(I,K,L)<0.) THEN 416 420 SM(I,K,L)=SM(I,K,L)+FM(I,K) 417 421 ALF(I,K)=FM(I,K)/SM(I,K,L) … … 427 431 ALF3(I,K)=ALF1(I,K)*ALF(I,K) 428 432 C 429 320 CONTINUE430 32 CONTINUE433 END DO 434 END DO 431 435 c print*,'ap ADVYP 32' 432 436 C 433 DO 33JV=1,NTRA434 DO 33K=1,LAT-1437 DO JV=1,NTRA 438 DO K=1,LAT-1 435 439 KP=K+1 436 DO 330I=1,LON437 C 438 IF(VGRI(I,K,L) .LT.0.) THEN440 DO I=1,LON 441 C 442 IF(VGRI(I,K,L)<0.) THEN 439 443 C 440 444 TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV) … … 474 478 ENDIF 475 479 C 476 330 CONTINUE 477 33 CONTINUE 480 END DO 481 END DO 482 END DO 478 483 c print*,'ap ADVYP 33' 479 484 C … … 483 488 C 484 489 SM0=0. 485 DO 40JV=1,NTRA490 DO JV=1,NTRA 486 491 S00(JV)=0. 487 40 CONTINUE488 C 489 DO 41I=1,LON490 C 491 IF(VGRI(I,K,L) .GE.0.) THEN492 END DO 493 C 494 DO I=1,LON 495 C 496 IF(VGRI(I,K,L)>=0.) THEN 492 497 FM(I,K)=VGRI(I,K,L)*DTY 493 498 ALF(I,K)=FM(I,K)/SM(I,K,L) … … 503 508 ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K) 504 509 C 505 41 CONTINUE510 END DO 506 511 c print*,'ap ADVYP 41' 507 512 C 508 DO 42JV=1,NTRA509 DO 420I=1,LON510 C 511 IF(VGRI(I,K,L) .GE.0.) THEN513 DO JV=1,NTRA 514 DO I=1,LON 515 C 516 IF(VGRI(I,K,L)>=0.) THEN 512 517 F0 (I,K,JV)=ALF(I,K)* ( S0(I,K,L,JV)+ALF1(I,K)* 513 518 + ( SY(I,K,L,JV)+ALF2(I,K)*SYY(I,K,L,JV) ) ) … … 527 532 ENDIF 528 533 C 529 420 CONTINUE530 42 CONTINUE534 END DO 535 END DO 531 536 c print*,'ap ADVYP 42' 532 537 C 533 DO 43I=1,LON534 IF(VGRI(I,K,L) .LT.0.) THEN538 DO I=1,LON 539 IF(VGRI(I,K,L)<0.) THEN 535 540 FM(I,K)=-VGRI(I,K,L)*DTY 536 541 ALF(I,K)=FM(I,K)/SM0 537 542 ENDIF 538 43 CONTINUE543 END DO 539 544 c print*,'ap ADVYP 43' 540 545 C 541 DO 44JV=1,NTRA542 DO 440I=1,LON543 IF(VGRI(I,K,L) .LT.0.) THEN546 DO JV=1,NTRA 547 DO I=1,LON 548 IF(VGRI(I,K,L)<0.) THEN 544 549 F0(I,K,JV)=ALF(I,K)*S00(JV) 545 550 ENDIF 546 440 CONTINUE547 44 CONTINUE551 END DO 552 END DO 548 553 C 549 554 C puts the temporary moments Fi into appropriate neighboring boxes 550 555 C 551 DO 45I=1,LON552 C 553 IF(VGRI(I,K,L) .LT.0.) THEN556 DO I=1,LON 557 C 558 IF(VGRI(I,K,L)<0.) THEN 554 559 SM(I,K,L)=SM(I,K,L)+FM(I,K) 555 560 ALF(I,K)=FM(I,K)/SM(I,K,L) … … 562 567 ALF3(I,K)=ALF1(I,K)*ALF(I,K) 563 568 C 564 45 CONTINUE569 END DO 565 570 c print*,'ap ADVYP 45' 566 571 C 567 DO 46JV=1,NTRA568 DO 460I=1,LON569 C 570 IF(VGRI(I,K,L) .LT.0.) THEN572 DO JV=1,NTRA 573 DO I=1,LON 574 C 575 IF(VGRI(I,K,L)<0.) THEN 571 576 C 572 577 TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV) … … 580 585 ENDIF 581 586 C 582 460 CONTINUE583 46 CONTINUE587 END DO 588 END DO 584 589 c print*,'ap ADVYP 46' 585 590 C 586 1 CONTINUE591 END DO 587 592 588 593 c-------------------------------------------------- -
LMDZ6/trunk/libf/dyn3d_common/advz.F
r4593 r5079 117 117 C Conversion du flux de masse en kg.s-1 118 118 119 DO 500l = 1,llm120 DO 500j = 1,jjp1121 DO 500i = 1,iip1119 DO l = 1,llm 120 DO j = 1,jjp1 121 DO i = 1,iip1 122 122 c wgri (i,j,llm+1-l) = w (i,j,l) / g 123 123 wgri (i,j,llm+1-l) = w (i,j,l) … … 125 125 c wgri (i,j,l) = 0.1 ! w (i,j,l) 126 126 c wgri (i,j,llm) = 0. ! a detruire ult. 127 500 CONTINUE 127 END DO 128 END DO 129 END DO 128 130 DO j = 1,jjp1 129 131 DO i = 1,iip1 … … 137 139 C boucle sur les latitudes 138 140 C 139 DO 1K=1,LAT141 DO K=1,LAT 140 142 C 141 143 C place limits on appropriate moments before transport … … 144 146 IF(.NOT.LIMIT) GO TO 101 145 147 C 146 DO 10JV=1,NTRA147 DO 10L=1,NIV148 DO 100I=1,LON148 DO JV=1,NTRA 149 DO L=1,NIV 150 DO I=1,LON 149 151 sz(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.), 150 152 + ABS(sz(I,K,L,JV))),sz(I,K,L,JV)) 151 100 CONTINUE 152 10 CONTINUE 153 END DO 154 END DO 155 END DO 153 156 C 154 157 101 CONTINUE … … 162 165 C 2- reajusts moments remaining in the box 163 166 C 164 DO 11L=1,NIV-1167 DO L=1,NIV-1 165 168 LP=L+1 166 169 C 167 DO 110I=1,LON168 C 169 IF(WGRI(I,K,L) .LT.0.) THEN170 DO I=1,LON 171 C 172 IF(WGRI(I,K,L)<0.) THEN 170 173 FM(I,L)=-WGRI(I,K,L)*DTZ 171 174 ALF(I)=FM(I,L)/SM(I,K,LP) … … 181 184 ALF1Q(I)=ALF1(I)*ALF1(I) 182 185 C 183 110 CONTINUE184 C 185 DO 111JV=1,NTRA186 DO 1110I=1,LON187 C 188 IF(WGRI(I,K,L) .LT.0.) THEN186 END DO 187 C 188 DO JV=1,NTRA 189 DO I=1,LON 190 C 191 IF(WGRI(I,K,L)<0.) THEN 189 192 C 190 193 F0(I,L,JV)=ALF (I)*( S0(I,K,LP,JV)-ALF1(I)*sz(I,K,LP,JV) ) … … 212 215 ENDIF 213 216 C 214 1110 CONTINUE215 111 CONTINUE216 C 217 11 CONTINUE217 END DO 218 END DO 219 C 220 END DO 218 221 C 219 222 C puts the temporary moments Fi into appropriate neighboring boxes 220 223 C 221 DO 12L=1,NIV-1224 DO L=1,NIV-1 222 225 LP=L+1 223 226 C 224 DO 120I=1,LON225 C 226 IF(WGRI(I,K,L) .LT.0.) THEN227 DO I=1,LON 228 C 229 IF(WGRI(I,K,L)<0.) THEN 227 230 SM(I,K,L)=SM(I,K,L)+FM(I,L) 228 231 ALF(I)=FM(I,L)/SM(I,K,L) … … 236 239 ALF1Q(I)=ALF1(I)*ALF1(I) 237 240 C 238 120 CONTINUE239 C 240 DO 121JV=1,NTRA241 DO 1210I=1,LON242 C 243 IF(WGRI(I,K,L) .LT.0.) THEN241 END DO 242 C 243 DO JV=1,NTRA 244 DO I=1,LON 245 C 246 IF(WGRI(I,K,L)<0.) THEN 244 247 C 245 248 TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV) … … 260 263 ENDIF 261 264 C 262 1210 CONTINUE263 121 CONTINUE264 C 265 12 CONTINUE265 END DO 266 END DO 267 C 268 END DO 266 269 C 267 270 C fin de la boucle principale sur les latitudes 268 271 C 269 1 CONTINUE272 END DO 270 273 C 271 274 C------------------------------------------------------------- -
LMDZ6/trunk/libf/dyn3d_common/advzp.F
r2600 r5079 135 135 C Conversion des flux de masses en kg 136 136 137 DO 500l = 1,llm138 DO 500j = 1,jjp1139 DO 500i = 1,iip1137 DO l = 1,llm 138 DO j = 1,jjp1 139 DO i = 1,iip1 140 140 wgri (i,j,llm+1-l) = w (i,j,l) 141 500 CONTINUE 141 END DO 142 END DO 143 END DO 142 144 do j=1,jjp1 143 145 do i=1,iip1 … … 154 156 C boucle sur les latitudes 155 157 C 156 DO 1K=1,LAT158 DO K=1,LAT 157 159 C 158 160 C place limits on appropriate moments before transport … … 161 163 IF(.NOT.LIMIT) GO TO 101 162 164 C 163 DO 10JV=1,NTRA164 DO 10L=1,NIV165 DO 100I=1,LON166 IF(S0(I,K,L,JV) .GT.0.) THEN165 DO JV=1,NTRA 166 DO L=1,NIV 167 DO I=1,LON 168 IF(S0(I,K,L,JV)>0.) THEN 167 169 SLPMAX=S0(I,K,L,JV) 168 170 S1MAX =1.5*SLPMAX … … 180 182 SYZ(I,K,L,JV)=0. 181 183 ENDIF 182 100 CONTINUE 183 10 CONTINUE 184 END DO 185 END DO 186 END DO 184 187 C 185 188 101 CONTINUE … … 193 196 C 2- reajusts moments remaining in the box 194 197 C 195 DO 11L=1,NIV-1198 DO L=1,NIV-1 196 199 LP=L+1 197 200 C 198 DO 110I=1,LON199 C 200 IF(WGRI(I,K,L) .LT.0.) THEN201 DO I=1,LON 202 C 203 IF(WGRI(I,K,L)<0.) THEN 201 204 FM(I,L)=-WGRI(I,K,L)*DTZ 202 205 ALF(I)=FM(I,L)/SM(I,K,LP) … … 215 218 ALF4 (I)=ALF1(I)*ALF1Q(I) 216 219 C 217 110 CONTINUE218 C 219 DO 111JV=1,NTRA220 DO 1110I=1,LON221 C 222 IF(WGRI(I,K,L) .LT.0.) THEN220 END DO 221 C 222 DO JV=1,NTRA 223 DO I=1,LON 224 C 225 IF(WGRI(I,K,L)<0.) THEN 223 226 C 224 227 F0 (I,L,JV)=ALF (I)* ( S0(I,K,LP,JV)-ALF1(I)* … … 273 276 ENDIF 274 277 C 275 1110 CONTINUE276 111 CONTINUE277 C 278 11 CONTINUE278 END DO 279 END DO 280 C 281 END DO 279 282 C 280 283 C puts the temporary moments Fi into appropriate neighboring boxes 281 284 C 282 DO 12L=1,NIV-1285 DO L=1,NIV-1 283 286 LP=L+1 284 287 C 285 DO 120I=1,LON286 C 287 IF(WGRI(I,K,L) .LT.0.) THEN288 DO I=1,LON 289 C 290 IF(WGRI(I,K,L)<0.) THEN 288 291 SM(I,K,L)=SM(I,K,L)+FM(I,L) 289 292 ALF(I)=FM(I,L)/SM(I,K,L) … … 299 302 ALF3(I)=ALF1(I)-ALF(I) 300 303 C 301 120 CONTINUE302 C 303 DO 121JV=1,NTRA304 DO 1210I=1,LON305 C 306 IF(WGRI(I,K,L) .LT.0.) THEN304 END DO 305 C 306 DO JV=1,NTRA 307 DO I=1,LON 308 C 309 IF(WGRI(I,K,L)<0.) THEN 307 310 C 308 311 TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV) … … 342 345 ENDIF 343 346 C 344 1210 CONTINUE345 121 CONTINUE346 C 347 12 CONTINUE347 END DO 348 END DO 349 C 350 END DO 348 351 C 349 352 C fin de la boucle principale sur les latitudes 350 353 C 351 1 CONTINUE354 END DO 352 355 C 353 356 DO l = 1,llm -
LMDZ6/trunk/libf/dyn3d_common/extrapol.F
r1952 r5079 58 58 200 CONTINUE 59 59 incre = incre + 1 60 DO 99999j = 1, kylat61 DO 99999i = 1, kxlon62 IF (pfild(i,j) .GT.zwmsk) THEN60 DO j = 1, kylat 61 DO i = 1, kxlon 62 IF (pfild(i,j)> zwmsk) THEN 63 63 pwork(i,j) = pfild(i,j) 64 64 inbor = 0 … … 89 89 C 90 90 C* Correct latitude bounds if southernmost or northernmost points 91 IF (j .EQ.1) ideb = 492 IF (j .EQ.kylat) ifin = 691 IF (j == 1) ideb = 4 92 IF (j == kylat) ifin = 6 93 93 C 94 94 C* Account for periodicity in longitude 95 95 C 96 96 IF (ldper) THEN 97 IF (i .EQ.kxlon) THEN97 IF (i == kxlon) THEN 98 98 ix(3) = 1 99 99 ix(6) = 1 100 100 ix(9) = 1 101 ELSE IF (i .EQ.1) THEN101 ELSE IF (i == 1) THEN 102 102 ix(1) = kxlon 103 103 ix(4) = kxlon … … 105 105 ENDIF 106 106 ELSE 107 IF (i .EQ.1) THEN107 IF (i == 1) THEN 108 108 ix(1) = i 109 109 ix(2) = i + 1 … … 113 113 ix(6) = i + 1 114 114 ENDIF 115 IF (i .EQ.kxlon) THEN115 IF (i == kxlon) THEN 116 116 ix(1) = i -1 117 117 ix(2) = i … … 122 122 ENDIF 123 123 C 124 IF (i .EQ. 1 .OR. i .EQ.kxlon) THEN124 IF (i == 1 .OR. i == kxlon) THEN 125 125 jy(1) = MAX (1,j-1) 126 126 jy(2) = MAX (1,j-1) … … 132 132 ideb = 1 133 133 ifin = 6 134 IF (j .EQ.1) ideb = 3135 IF (j .EQ.kylat) ifin = 4134 IF (j == 1) ideb = 3 135 IF (j == kylat) ifin = 4 136 136 ENDIF 137 137 ENDIF ! end for ldper test … … 139 139 C* Find unmasked neighbors 140 140 C 141 DO 230k = ideb, ifin141 DO k = ideb, ifin 142 142 zmask(k) = 0. 143 143 ilon = ix(k) 144 144 jlat = jy(k) 145 IF (pfild(ilon,jlat) .LT.zwmsk) THEN145 IF (pfild(ilon,jlat) < zwmsk) THEN 146 146 zmask(k) = 1. 147 147 inbor = inbor + 1 148 148 ENDIF 149 230 CONTINUE149 END DO 150 150 C 151 151 C* Not enough points around point P are unmasked; interpolation on P 152 152 C will be done in a future call to extrap. 153 153 C 154 IF (inbor .GE.knbor) THEN154 IF (inbor >= knbor) THEN 155 155 pwork(i,j) = 0. 156 156 DO k = ideb, ifin … … 163 163 C 164 164 ENDIF 165 99999 CONTINUE 165 END DO 166 END DO 166 167 C 167 168 C* 3. Writing back unmasked field in pfild … … 176 177 DO j = 1, kylat 177 178 DO i = 1, kxlon 178 IF (pwork(i,j) .GT.zwmsk) idoit = idoit + 1179 IF (pwork(i,j) > zwmsk) idoit = idoit + 1 179 180 pfild(i,j) = pwork(i,j) 180 181 ENDDO 181 182 ENDDO 182 183 c 183 IF (idoit .ne.0) GOTO 200184 IF (idoit /= 0) GOTO 200 184 185 ccc PRINT*, "Number of extrapolation steps incre =", incre 185 186 c -
LMDZ6/trunk/libf/dyn3d_common/ppm3d.F
r2197 r5079 68 68 implicit none 69 69 70 c rajout de d éclarations70 c rajout de d�clarations 71 71 c integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd 72 72 c integer iu,iiu,j2,jmr,js0,jt … … 315 315 C 316 316 C *********** Initialization ********************** 317 if(NSTEP .eq.1) then317 if(NSTEP==1) then 318 318 c 319 319 write(6,*) '------------------------------------ ' … … 325 325 C 326 326 C controles sur les parametres 327 if(NLAY .LT.6) then327 if(NLAY<6) then 328 328 write(6,*) 'NLAY must be >= 6' 329 329 stop 330 330 endif 331 if (JNP .LT.NLAY) then331 if (JNP<NLAY) then 332 332 write(6,*) 'JNP must be >= NLAY' 333 333 stop 334 334 endif 335 335 IMRD2=mod(IMR,2) 336 if (j1 .eq.2.and.IMRD2.NE.0) then336 if (j1==2.and.IMRD2/=0) then 337 337 write(6,*) 'if j1=2 IMR must be an even integer' 338 338 stop … … 340 340 341 341 C 342 if(Jmax .lt.JNP .or. Kmax.lt.NLAY) then342 if(Jmax<JNP .or. Kmax<NLAY) then 343 343 write(6,*) 'Jmax or Kmax is too small' 344 344 stop … … 354 354 DP = PI / REAL(JMR) 355 355 C 356 if(IGD .eq.0) then356 if(IGD==0) then 357 357 C Compute analytic cosine at cell edges 358 358 call cosa(cosp,cose,JNP,PI,DP) … … 362 362 endif 363 363 C 364 do 15 J=2,JMR 365 15 acosp(j) = 1. / cosp(j) 364 do J=2,JMR 365 acosp(j) = 1. / cosp(j) 366 END DO 366 367 C 367 368 C Inverse of the Scaled polar cap area. … … 372 373 endif 373 374 C 374 if(NDT0 .ne.NDT) then375 if(NDT0 /= NDT) then 375 376 DT = NDT 376 377 NDT0 = NDT 377 378 378 if(Umax .lt.180.) then379 if(Umax < 180.) then 379 380 write(6,*) 'Umax may be too small!' 380 381 endif … … 382 383 MaxDT = DP*AE / abs(Umax) + 0.5 383 384 write(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT 384 if(MaxDT .lt.abs(NDT)) then385 if(MaxDT < abs(NDT)) then 385 386 write(6,*) 'Warning!!! NDT maybe too large!' 386 387 endif 387 388 C 388 if(CR1 .ge.0.95) then389 if(CR1>=0.95) then 389 390 JS0 = 0 390 391 JN0 = 0 … … 429 430 430 431 C 431 if(j1 .ne.2) then432 DO 40IC=1,NC433 DO 40L=1,NLAY434 DO 40I=1,IMR432 if(j1/=2) then 433 DO IC=1,NC 434 DO L=1,NLAY 435 DO I=1,IMR 435 436 Q(I, 2,L,IC) = Q(I, 1,L,IC) 436 40 Q(I,JMR,L,IC) = Q(I,JNP,L,IC) 437 Q(I,JMR,L,IC) = Q(I,JNP,L,IC) 438 END DO 439 END DO 440 END DO 437 441 endif 438 442 C 439 443 C Compute "tracer density" 440 DO 550 IC=1,NC 441 DO 44 k=1,NLAY 442 DO 44 j=1,JNP 443 DO 44 i=1,IMR 444 44 DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k) 445 550 continue 446 C 447 do 1500 k=1,NLAY 448 C 449 if(IGD.eq.0) then 444 DO IC=1,NC 445 DO k=1,NLAY 446 DO j=1,JNP 447 DO i=1,IMR 448 DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k) 449 END DO 450 END DO 451 END DO 452 END DO 453 C 454 do k=1,NLAY 455 C 456 if(IGD==0) then 450 457 C Convert winds on A-Grid to Courant # on C-Grid. 451 458 call A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5) 452 459 else 453 460 C Convert winds on C-grid to Courant # 454 do 45 j=j1,j2 455 do 45 i=2,IMR 456 45 CRX(i,J) = dtdx(j)*U(i-1,j,k) 461 do j=j1,j2 462 do i=2,IMR 463 CRX(i,J) = dtdx(j)*U(i-1,j,k) 464 END DO 465 END DO 457 466 458 467 C 459 do 50 j=j1,j2 460 50 CRX(1,J) = dtdx(j)*U(IMR,j,k) 461 C 462 do 55 i=1,IMR*JMR 463 55 CRY(i,2) = DTDY*V(i,1,k) 468 do j=j1,j2 469 CRX(1,J) = dtdx(j)*U(IMR,j,k) 470 END DO 471 C 472 do i=1,IMR*JMR 473 CRY(i,2) = DTDY*V(i,1,k) 474 END DO 464 475 endif 465 476 C … … 470 481 do j=JS0,j1+1,-1 471 482 do i=1,IMR 472 if(abs(CRX(i,j)) .GT.1.) then483 if(abs(CRX(i,j))>1.) then 473 484 JS = j 474 485 go to 2222 … … 480 491 do j=JN0,j2-1 481 492 do i=1,IMR 482 if(abs(CRX(i,j)) .GT.1.) then493 if(abs(CRX(i,j))>1.) then 483 494 JN = j 484 495 go to 2233 … … 488 499 2233 continue 489 500 C 490 if(j1 .ne.2) then ! Enlarged polar cap.501 if(j1/=2) then ! Enlarged polar cap. 491 502 do i=1,IMR 492 503 DPI(i, 2,k) = 0. … … 505 516 enddo 506 517 C 507 do 95 j=j1,j2 508 DO 95 i=1,IMR 509 95 DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j) 518 do j=j1,j2 519 DO i=1,IMR 520 DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j) 521 END DO 522 END DO 510 523 C 511 524 C Poles … … 536 549 enddo 537 550 C 538 do 110 j=j1,j2 539 DO 110 i=1,IMR 540 110 xmass(i,j) = PU(i,j)*CRX(i,j) 541 C 542 DO 120 j=j1,j2 543 DO 120 i=1,IMR-1 544 120 DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j) 545 C 546 DO 130 j=j1,j2 547 130 DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j) 551 do j=j1,j2 552 DO i=1,IMR 553 xmass(i,j) = PU(i,j)*CRX(i,j) 554 END DO 555 END DO 556 C 557 DO j=j1,j2 558 DO i=1,IMR-1 559 DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j) 560 END DO 561 END DO 562 C 563 DO j=j1,j2 564 DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j) 565 END DO 548 566 C 549 567 DO j=j1,j2 … … 569 587 enddo 570 588 C 571 if(j1 .eq.2) then589 if(j1==2) then 572 590 IMH = IMR/2 573 591 do i=1,IMH … … 582 600 C 583 601 C ****6***0*********0*********0*********0*********0*********0**********72 584 do 1000IC=1,NC602 do IC=1,NC 585 603 C 586 604 do i=1,IMJM … … 590 608 C 591 609 C E-W advective cross term 592 do 250j=J1,J2593 if(J .GT.JS .and. J.LT.JN) GO TO 250610 do j=J1,J2 611 if(J>JS .and. J<JN) GO TO 250 594 612 C 595 613 do i=1,IMR … … 602 620 enddo 603 621 C 604 DO 230i=1,IMR622 DO i=1,IMR 605 623 iu = UA(i,j) 606 624 ru = UA(i,j) - iu 607 625 iiu = i-iu 608 if(UA(i,j) .GE.0.) then626 if(UA(i,j)>=0.) then 609 627 wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu)) 610 628 else … … 612 630 endif 613 631 wk1(i,j,1) = wk1(i,j,1) - qtmp(i) 614 230 continue 632 END DO 615 633 250 continue 616 C 617 if(JN.ne.0) then 634 END DO 635 C 636 if(JN/=0) then 618 637 do j=JS+1,JN-1 619 638 C … … 645 664 if(cross) then 646 665 C Add cross terms in the vertical direction. 647 if(IORD .GE.2) then666 if(IORD >= 2) then 648 667 iad = 2 649 668 else … … 651 670 endif 652 671 C 653 if(JORD .GE.2) then672 if(JORD >= 2) then 654 673 jad = 2 655 674 else … … 671 690 & DC2,ymass,WK1(1,1,3),wk1(1,1,4),WK1(1,1,5),WK1(1,1,6),JORD) 672 691 C 673 1000 continue 674 1500 continue 692 END DO 693 END DO 675 694 C 676 695 C ******* Compute vertical mass flux (same unit as PS) *********** … … 678 697 C 1st step: compute total column mass CONVERGENCE. 679 698 C 680 do 320 j=1,JNP 681 do 320 i=1,IMR 682 320 CRY(i,j) = DPI(i,j,1) 683 C 684 do 330 k=2,NLAY 685 do 330 j=1,JNP 686 do 330 i=1,IMR 699 do j=1,JNP 700 do i=1,IMR 701 CRY(i,j) = DPI(i,j,1) 702 END DO 703 END DO 704 C 705 do k=2,NLAY 706 do j=1,JNP 707 do i=1,IMR 687 708 CRY(i,j) = CRY(i,j) + DPI(i,j,k) 688 330 continue 689 C 690 do 360 j=1,JNP 691 do 360 i=1,IMR 709 END DO 710 END DO 711 END DO 712 C 713 do j=1,JNP 714 do i=1,IMR 692 715 C 693 716 C 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption. … … 700 723 W(i,j,1) = DPI(i,j,1) - DBK(1)*CRY(i,j) 701 724 W(i,j,NLAY) = 0. 702 360 continue 703 C 704 do 370 k=2,NLAY-1 705 do 370 j=1,JNP 706 do 370 i=1,IMR 725 END DO 726 END DO 727 C 728 do k=2,NLAY-1 729 do j=1,JNP 730 do i=1,IMR 707 731 W(i,j,k) = W(i,j,k-1) + DPI(i,j,k) - DBK(k)*CRY(i,j) 708 370 continue 709 C 710 DO 380 k=1,NLAY 711 DO 380 j=1,JNP 712 DO 380 i=1,IMR 732 END DO 733 END DO 734 END DO 735 C 736 DO k=1,NLAY 737 DO j=1,JNP 738 DO i=1,IMR 713 739 delp2(i,j,k) = DAP(k) + DBK(k)*PS2(i,j) 714 380 continue 740 END DO 741 END DO 742 END DO 715 743 C 716 744 KRD = max(3, KORD) 717 do 4000IC=1,NC745 do IC=1,NC 718 746 C 719 747 C****6***0*********0*********0*********0*********0*********0**********72 … … 738 766 enddo 739 767 C 740 if(j1 .ne.2) then741 DO 400k=1,NLAY742 DO 400I=1,IMR743 c j=1 c'est le p ôle Sud, j=JNP c'est le pôle Nord768 if(j1/=2) then 769 DO k=1,NLAY 770 DO I=1,IMR 771 c j=1 c'est le p�le Sud, j=JNP c'est le p�le Nord 744 772 Q(I, 2,k,IC) = Q(I, 1,k,IC) 745 773 Q(I,JMR,k,IC) = Q(I,JNP,k,IC) 746 400 CONTINUE 747 endif 748 4000 continue 749 C 750 if(j1.ne.2) then 751 DO 5000 k=1,NLAY 752 DO 5000 i=1,IMR 774 END DO 775 END DO 776 endif 777 END DO 778 C 779 if(j1/=2) then 780 DO k=1,NLAY 781 DO i=1,IMR 753 782 W(i, 2,k) = W(i, 1,k) 754 783 W(i,JMR,k) = W(i,JNP,k) 755 5000 continue 784 END DO 785 END DO 756 786 endif 757 787 C … … 785 815 C ****6***0*********0*********0*********0*********0*********0**********72 786 816 C 787 do 1000k=1,NLAYM1788 do 1000i=1,IMJM817 do k=1,NLAYM1 818 do i=1,IMJM 789 819 DQDT(i,1,k) = P(i,1,k+1) - P(i,1,k) 790 1000 continue 791 C 792 DO 1220 k=2,NLAYM1 793 DO 1220 I=1,IMJM 820 END DO 821 END DO 822 C 823 DO k=2,NLAYM1 824 DO I=1,IMJM 794 825 c0 = delp(i,1,k) / (delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1)) 795 826 c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k)) … … 799 830 Qmin = P(i,1,k) - min(P(i,1,k-1),P(i,1,k),P(i,1,k+1)) 800 831 DC(i,1,k) = sign(min(abs(tmp),Qmax,Qmin), tmp) 801 1220 CONTINUE 832 END DO 833 END DO 802 834 803 835 C … … 806 838 C ****6***0*********0*********0*********0*********0*********0**********72 807 839 C 808 DO 2000j=1,JNP809 if((j .eq.2 .or. j.eq.JMR) .and. j1.ne.2) goto 2000840 DO j=1,JNP 841 if((j==2 .or. j==JMR) .and. j1/=2) goto 2000 810 842 C 811 843 DO k=1,NLAY … … 828 860 C 829 861 C First guess top edge value 830 DO 10i=1,IMR862 DO i=1,IMR 831 863 C three-cell PPM 832 864 C Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp … … 840 872 C 841 873 C Check if change sign 842 if(wk1(i,1)*AL(i,1) .le.0.) then874 if(wk1(i,1)*AL(i,1)<=0.) then 843 875 AL(i,1) = 0. 844 876 flux(i,1) = 0. … … 846 878 flux(i,1) = wk1(i,1) - AL(i,1) 847 879 endif 848 10 continue 880 END DO 849 881 C 850 882 C Bottom 851 DO 15i=1,IMR883 DO i=1,IMR 852 884 C 2-cell PPM with zero gradient right at the surface 853 885 C … … 856 888 AR(i,NLAY) = wk1(i,NLAY) + fct 857 889 AL(i,NLAY) = wk1(i,NLAY) - (fct+fct) 858 if(wk1(i,NLAY)*AR(i,NLAY) .le.0.) AR(i,NLAY) = 0.890 if(wk1(i,NLAY)*AR(i,NLAY)<=0.) AR(i,NLAY) = 0. 859 891 flux(i,NLAY) = AR(i,NLAY) - wk1(i,NLAY) 860 15 continue 892 END DO 861 893 862 894 C … … 865 897 C****6***0*********0*********0*********0*********0*********0**********72 866 898 C 867 DO 14k=3,NLAYM1868 DO 12i=1,IMR899 DO k=3,NLAYM1 900 DO i=1,IMR 869 901 c1 = DQDT(i,j,k-1)*wk2(i,k-1) / (wk2(i,k-1)+wk2(i,k)) 870 902 c2 = 2. / (wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1)) … … 875 907 & wk2(i,k-1)*A1*flux(i,k) ) 876 908 C print *,'AL1',i,k, AL(i,k) 877 12 CONTINUE 878 14 continue 879 C 880 do 20i=1,IMR*NLAYM1909 END DO 910 END DO 911 C 912 do i=1,IMR*NLAYM1 881 913 AR(i,1) = AL(i,2) 882 914 C print *,'AR1',i,AR(i,1) 883 20 continue 884 C 885 do 30i=1,IMR*NLAY915 END DO 916 C 917 do i=1,IMR*NLAY 886 918 A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1))) 887 919 C print *,'A61',i,A6(i,1) 888 30 continue 920 END DO 889 921 C 890 922 C****6***0*********0*********0*********0*********0*********0**********72 … … 895 927 C 896 928 C Interior depending on KORD 897 if(LMT .LE.2)929 if(LMT<=2) 898 930 & call lmtppm(flux(1,2),A6(1,2),AR(1,2),AL(1,2),wk1(1,2), 899 931 & IMR*(NLAY-2),LMT) … … 901 933 C****6***0*********0*********0*********0*********0*********0**********72 902 934 C 903 DO 140i=1,IMR*NLAYM1904 IF(wz2(i,1) .GT.0.) then935 DO i=1,IMR*NLAYM1 936 IF(wz2(i,1)>0.) then 905 937 CM = wz2(i,1) / wk2(i,1) 906 938 flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM)) … … 912 944 C print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23 913 945 endif 914 140 continue 915 C 916 DO 250i=1,IMR*NLAYM1946 END DO 947 C 948 DO i=1,IMR*NLAYM1 917 949 flux(i,2) = wz2(i,1) * flux(i,2) 918 250 continue 919 C 920 do 350i=1,IMR950 END DO 951 C 952 do i=1,IMR 921 953 DQ(i,j, 1) = DQ(i,j, 1) - flux(i, 2) 922 954 DQ(i,j,NLAY) = DQ(i,j,NLAY) + flux(i,NLAY) 923 350 continue 924 C 925 do 360 k=2,NLAYM1 926 do 360 i=1,IMR 927 360 DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1) 955 END DO 956 C 957 do k=2,NLAYM1 958 do i=1,IMR 959 DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1) 960 END DO 961 END DO 928 962 2000 continue 963 END DO 929 964 return 930 965 end … … 950 985 j2vl = j2-jvan 951 986 C 952 do 1310j=j1,j2987 do j=j1,j2 953 988 C 954 989 do i=1,IMR … … 956 991 enddo 957 992 C 958 if(j .ge.JN .or. j.le.JS) goto 2222993 if(j>=JN .or. j<=JS) goto 2222 959 994 C ************* Eulerian ********** 960 995 C … … 964 999 qtmp(IMP+1) = q(2,J) 965 1000 C 966 IF(IORD .eq.1 .or. j.eq.j1. or. j.eq.j2) THEN967 DO 1406i=1,IMR1001 IF(IORD==1 .or. j==j1 .or. j==j2) THEN 1002 DO i=1,IMR 968 1003 iu = REAL(i) - uc(i,j) 969 1406 fx1(i) = qtmp(iu) 1004 fx1(i) = qtmp(iu) 1005 END DO 970 1006 ELSE 971 1007 call xmist(IMR,IML,Qtmp,DC) 972 1008 DC(0) = DC(IMR) 973 1009 C 974 if(IORD .eq.2 .or. j.le.j1vl .or. j.ge.j2vl) then975 DO 1408i=1,IMR1010 if(IORD==2 .or. j<=j1vl .or. j>=j2vl) then 1011 DO i=1,IMR 976 1012 iu = REAL(i) - uc(i,j) 977 1408 fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j)) 1013 fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j)) 1014 END DO 978 1015 else 979 1016 call fxppm(IMR,IML,UC(1,j),Qtmp,DC,fx1,IORD) … … 982 1019 ENDIF 983 1020 C 984 DO 1506 i=1,IMR 985 1506 fx1(i) = fx1(i)*xmass(i,j) 1021 DO i=1,IMR 1022 fx1(i) = fx1(i)*xmass(i,j) 1023 END DO 986 1024 C 987 1025 goto 1309 … … 996 1034 enddo 997 1035 C 998 IF(IORD .eq.1 .or. j.eq.j1. or. j.eq.j2) THEN999 DO 1306i=1,IMR1036 IF(IORD==1 .or. j==j1 .or. j==j2) THEN 1037 DO i=1,IMR 1000 1038 itmp = INT(uc(i,j)) 1001 1039 ISAVE(i) = i - itmp 1002 1040 iu = i - uc(i,j) 1003 1306 fx1(i) = (uc(i,j) - itmp)*qtmp(iu) 1041 fx1(i) = (uc(i,j) - itmp)*qtmp(iu) 1042 END DO 1004 1043 ELSE 1005 1044 call xmist(IMR,IML,Qtmp,DC) … … 1010 1049 enddo 1011 1050 C 1012 DO 1307i=1,IMR1051 DO i=1,IMR 1013 1052 itmp = INT(uc(i,j)) 1014 1053 rut = uc(i,j) - itmp 1015 1054 ISAVE(i) = i - itmp 1016 1055 iu = i - uc(i,j) 1017 1307 fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut)) 1056 fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut)) 1057 END DO 1018 1058 ENDIF 1019 1059 C 1020 do 1308i=1,IMR1021 IF(uc(i,j) .GT.1.) then1060 do i=1,IMR 1061 IF(uc(i,j)>1.) then 1022 1062 CDIR$ NOVECTOR 1023 1063 do ist = ISAVE(i),i-1 1024 1064 fx1(i) = fx1(i) + qtmp(ist) 1025 1065 enddo 1026 elseIF(uc(i,j) .LT.-1.) then1066 elseIF(uc(i,j)<-1.) then 1027 1067 do ist = i,ISAVE(i)-1 1028 1068 fx1(i) = fx1(i) - qtmp(ist) … … 1030 1070 CDIR$ VECTOR 1031 1071 endif 1032 1308 continue 1072 END DO 1033 1073 do i=1,IMR 1034 1074 fx1(i) = PU(i,j)*fx1(i) … … 1038 1078 C 1039 1079 1309 fx1(IMP) = fx1(1) 1040 DO 1215 i=1,IMR 1041 1215 DQ(i,j) = DQ(i,j) + fx1(i)-fx1(i+1) 1080 DO i=1,IMR 1081 DQ(i,j) = DQ(i,j) + fx1(i)-fx1(i+1) 1082 END DO 1042 1083 C 1043 1084 C *************************************** 1044 1085 C 1045 1310 continue 1086 END DO 1046 1087 return 1047 1088 end … … 1079 1120 C endif 1080 1121 C 1081 DO 10 i=1,IMR 1082 10 AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3 1083 C 1084 do 20 i=1,IMR-1 1085 20 AR(i) = AL(i+1) 1122 DO i=1,IMR 1123 AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3 1124 END DO 1125 C 1126 do i=1,IMR-1 1127 AR(i) = AL(i+1) 1128 END DO 1086 1129 AR(IMR) = AL(1) 1087 1130 C 1088 do 30 i=1,IMR 1089 30 A6(i) = 3.*(p(i)+p(i) - (AL(i)+AR(i))) 1090 C 1091 if(LMT.LE.2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT) 1131 do i=1,IMR 1132 A6(i) = 3.*(p(i)+p(i) - (AL(i)+AR(i))) 1133 END DO 1134 C 1135 if(LMT<=2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT) 1092 1136 C 1093 1137 AL(0) = AL(IMR) … … 1096 1140 C 1097 1141 DO i=1,IMR 1098 IF(UT(i) .GT.0.) then1142 IF(UT(i)>0.) then 1099 1143 flux(i) = AR(i-1) + 0.5*UT(i)*(AL(i-1) - AR(i-1) + 1100 1144 & A6(i-1)*(1.-R23*UT(i)) ) … … 1115 1159 real :: tmp,pmax,pmin 1116 1160 C 1117 do 10i=1,IMR1161 do i=1,IMR 1118 1162 tmp = R24*(8.*(p(i+1) - p(i-1)) + p(i-2) - p(i+2)) 1119 1163 Pmax = max(P(i-1), p(i), p(i+1)) - p(i) 1120 1164 Pmin = p(i) - min(P(i-1), p(i), p(i+1)) 1121 10 DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp) 1165 DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp) 1166 END DO 1122 1167 return 1123 1168 end … … 1138 1183 len = IMR*(J2-J1+2) 1139 1184 C 1140 if(JORD .eq.1) then1141 DO 1000i=1,len1185 if(JORD==1) then 1186 DO i=1,len 1142 1187 JT = REAL(J1) - VC(i,J1) 1143 1000 fx(i,j1) = p(i,JT) 1188 fx(i,j1) = p(i,JT) 1189 END DO 1144 1190 else 1145 1191 1146 1192 call ymist(IMR,JNP,j1,P,DC2,4) 1147 1193 C 1148 if(JORD .LE.0 .or. JORD.GE.3) then1194 if(JORD<=0 .or. JORD>=3) then 1149 1195 1150 1196 call fyppm(VC,P,DC2,fx,IMR,JNP,j1,j2,A6,AR,AL,JORD) 1151 1197 1152 1198 else 1153 DO 1200i=1,len1199 DO i=1,len 1154 1200 JT = REAL(J1) - VC(i,J1) 1155 1200 fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT) 1156 endif 1157 endif 1158 C 1159 DO 1300 i=1,len 1160 1300 fx(i,j1) = fx(i,j1)*ymass(i,j1) 1161 C 1162 DO 1400 j=j1,j2 1163 DO 1400 i=1,IMR 1164 1400 DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j) 1201 fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT) 1202 END DO 1203 endif 1204 endif 1205 C 1206 DO i=1,len 1207 fx(i,j1) = fx(i,j1)*ymass(i,j1) 1208 END DO 1209 C 1210 DO j=j1,j2 1211 DO i=1,IMR 1212 DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j) 1213 END DO 1214 END DO 1165 1215 C 1166 1216 C Poles … … 1179 1229 enddo 1180 1230 C 1181 if(j1 .ne.2) then1231 if(j1/=2) then 1182 1232 do i=1,IMR 1183 1233 DQ(i, 2) = sum1 … … 1201 1251 IJM3 = IMR*(JMR-3) 1202 1252 C 1203 IF(ID .EQ.2) THEN1204 do 10i=1,IMR*(JMR-1)1253 IF(ID==2) THEN 1254 do i=1,IMR*(JMR-1) 1205 1255 tmp = 0.25*(p(i,3) - p(i,1)) 1206 1256 Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2) 1207 1257 Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3)) 1208 1258 DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp) 1209 10 CONTINUE 1259 END DO 1210 1260 ELSE 1211 do 12i=1,IMH1261 do i=1,IMH 1212 1262 C J=2 1213 1263 tmp = (8.*(p(i,3) - p(i,1)) + p(i+IMH,2) - p(i,4))*R24 … … 1220 1270 Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP)) 1221 1271 DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp) 1222 12 CONTINUE 1223 do 14i=IMH+1,IMR1272 END DO 1273 do i=IMH+1,IMR 1224 1274 C J=2 1225 1275 tmp = (8.*(p(i,3) - p(i,1)) + p(i-IMH,2) - p(i,4))*R24 … … 1232 1282 Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP)) 1233 1283 DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp) 1234 14 CONTINUE 1235 C 1236 do 15i=1,IJM31284 END DO 1285 C 1286 do i=1,IJM3 1237 1287 tmp = (8.*(p(i,4) - p(i,2)) + p(i,1) - p(i,5))*R24 1238 1288 Pmax = max(p(i,2),p(i,3),p(i,4)) - p(i,3) 1239 1289 Pmin = p(i,3) - min(p(i,2),p(i,3),p(i,4)) 1240 1290 DC(i,3) = sign(min(abs(tmp),Pmin,Pmax),tmp) 1241 15 CONTINUE 1291 END DO 1242 1292 ENDIF 1243 1293 C 1244 if(j1 .ne.2) then1294 if(j1/=2) then 1245 1295 do i=1,IMR 1246 1296 DC(i,1) = 0. … … 1250 1300 C Determine slopes in polar caps for scalars! 1251 1301 C 1252 do 13i=1,IMH1302 do i=1,IMH 1253 1303 C South 1254 1304 tmp = 0.25*(p(i,2) - p(i+imh,2)) … … 1261 1311 Pmin = p(i,JNP) - min(p(i+imh,JMR),p(i,jnp), p(i,JMR)) 1262 1312 DC(i,JNP) = sign(min(abs(tmp),Pmax,pmin),tmp) 1263 13 continue 1264 C 1265 do 25i=imh+1,IMR1313 END DO 1314 C 1315 do i=imh+1,IMR 1266 1316 DC(i, 1) = - DC(i-imh, 1) 1267 1317 DC(i,JNP) = - DC(i-imh,JNP) 1268 25 continue 1318 END DO 1269 1319 endif 1270 1320 return … … 1308 1358 LMT = JORD - 3 1309 1359 C 1310 DO 10 i=1,IMR*JMR1360 DO i=1,IMR*JMR 1311 1361 AL(i,2) = 0.5*(p(i,1)+p(i,2)) + (DC(i,1) - DC(i,2))*R3 1312 1362 AR(i,1) = AL(i,2) 1313 10 CONTINUE 1363 END DO 1314 1364 C 1315 1365 CPoles: … … 1331 1381 1332 1382 1333 do 30 i=1,len 1334 30 A6(i,j11) = 3.*(p(i,j11)+p(i,j11) - (AL(i,j11)+AR(i,j11))) 1335 C 1336 if(LMT.le.2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11) 1383 do i=1,len 1384 A6(i,j11) = 3.*(p(i,j11)+p(i,j11) - (AL(i,j11)+AR(i,j11))) 1385 END DO 1386 C 1387 if(LMT<=2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11) 1337 1388 & ,AL(1,j11),P(1,j11),len,LMT) 1338 1389 C 1339 1390 1340 DO 140i=1,IMJM11341 IF(VC(i,j1) .GT.0.) then1391 DO i=1,IMJM1 1392 IF(VC(i,j1)>0.) then 1342 1393 flux(i,j1) = AR(i,j11) + 0.5*VC(i,j1)*(AL(i,j11) - AR(i,j11) + 1343 1394 & A6(i,j11)*(1.-R23*VC(i,j1)) ) … … 1346 1397 & A6(i,j1)*(1.+R23*VC(i,j1))) 1347 1398 endif 1348 140 continue 1399 END DO 1349 1400 return 1350 1401 end … … 1378 1429 c write(*,*) 'toto 1' 1379 1430 C -------------------------------- 1380 IF(IAD .eq.2) then1431 IF(IAD==2) then 1381 1432 do j=j1-1,j2+1 1382 1433 do i=1,IMR … … 1395 1446 c write(*,*) 'toto 2' 1396 1447 C 1397 ELSEIF(IAD .eq.1) then1448 ELSEIF(IAD==1) then 1398 1449 do j=j1-1,j2+1 1399 1450 do i=1,imr … … 1404 1455 ENDIF 1405 1456 C 1406 if(j1 .ne.2) then1457 if(j1/=2) then 1407 1458 sum1 = 0. 1408 1459 sum2 = 0. … … 1448 1499 C 1449 1500 JMR = JNP-1 1450 do 1309j=j1,j21451 if(J .GT.JS .and. J.LT.JN) GO TO 13091501 do j=j1,j2 1502 if(J>JS .and. J<JN) GO TO 1309 1452 1503 C 1453 1504 do i=1,IMR … … 1460 1511 enddo 1461 1512 C 1462 IF(IAD .eq.2) THEN1513 IF(IAD==2) THEN 1463 1514 DO i=1,IMR 1464 1515 IP = NINT(UA(i,j)) … … 1469 1520 adx(i,j) = qtmp(ip) + ru*(a1*ru + b1) 1470 1521 enddo 1471 ELSEIF(IAD .eq.1) then1522 ELSEIF(IAD==1) then 1472 1523 DO i=1,IMR 1473 1524 iu = UA(i,j) 1474 1525 ru = UA(i,j) - iu 1475 1526 iiu = i-iu 1476 if(UA(i,j) .GE.0.) then1527 if(UA(i,j)>=0.) then 1477 1528 adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu)) 1478 1529 else … … 1486 1537 enddo 1487 1538 1309 continue 1539 END DO 1488 1540 C 1489 1541 C Eulerian upwind … … 1498 1550 qtmp(IMR+1) = p(1,J) 1499 1551 C 1500 IF(IAD .eq.2) THEN1552 IF(IAD==2) THEN 1501 1553 qtmp(-1) = p(IMR-1,J) 1502 1554 qtmp(IMR+2) = p(2,J) … … 1509 1561 adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1) 1510 1562 enddo 1511 ELSEIF(IAD .eq.1) then1563 ELSEIF(IAD==1) then 1512 1564 C 1st order 1513 1565 DO i=1,IMR … … 1518 1570 enddo 1519 1571 C 1520 if(j1 .ne.2) then1572 if(j1/=2) then 1521 1573 do i=1,IMR 1522 1574 adx(i, 2) = 0. … … 1554 1606 REAL da1,da2,a6da,fmin 1555 1607 C 1556 if(LMT .eq.0) then1608 if(LMT==0) then 1557 1609 C Full constraint 1558 do 100i=1,IM1559 if(DC(i) .eq.0.) then1610 do i=1,IM 1611 if(DC(i)==0.) then 1560 1612 AR(i) = p(i) 1561 1613 AL(i) = p(i) … … 1565 1617 da2 = da1**2 1566 1618 A6DA = A6(i)*da1 1567 if(A6DA .lt.-da2) then1619 if(A6DA < -da2) then 1568 1620 A6(i) = 3.*(AL(i)-p(i)) 1569 1621 AR(i) = AL(i) - A6(i) 1570 elseif(A6DA .gt.da2) then1622 elseif(A6DA > da2) then 1571 1623 A6(i) = 3.*(AR(i)-p(i)) 1572 1624 AL(i) = AR(i) - A6(i) 1573 1625 endif 1574 1626 endif 1575 100 continue 1576 elseif(LMT .eq.1) then1627 END DO 1628 elseif(LMT==1) then 1577 1629 C Semi-monotonic constraint 1578 do 150i=1,IM1579 if(abs(AR(i)-AL(i)) .GE.-A6(i)) go to 1501580 if(p(i) .lt.AR(i) .and. p(i).lt.AL(i)) then1630 do i=1,IM 1631 if(abs(AR(i)-AL(i)) >= -A6(i)) go to 150 1632 if(p(i)<AR(i) .and. p(i)<AL(i)) then 1581 1633 AR(i) = p(i) 1582 1634 AL(i) = p(i) 1583 1635 A6(i) = 0. 1584 elseif(AR(i) .gt.AL(i)) then1636 elseif(AR(i) > AL(i)) then 1585 1637 A6(i) = 3.*(AL(i)-p(i)) 1586 1638 AR(i) = AL(i) - A6(i) … … 1590 1642 endif 1591 1643 150 continue 1592 elseif(LMT.eq.2) then 1593 do 250 i=1,IM 1594 if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 250 1644 END DO 1645 elseif(LMT==2) then 1646 do i=1,IM 1647 if(abs(AR(i)-AL(i)) >= -A6(i)) go to 250 1595 1648 fmin = p(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12 1596 if(fmin .ge.0.) go to 2501597 if(p(i) .lt.AR(i) .and. p(i).lt.AL(i)) then1649 if(fmin>=0.) go to 250 1650 if(p(i)<AR(i) .and. p(i)<AL(i)) then 1598 1651 AR(i) = p(i) 1599 1652 AL(i) = p(i) 1600 1653 A6(i) = 0. 1601 elseif(AR(i) .gt.AL(i)) then1654 elseif(AR(i) > AL(i)) then 1602 1655 A6(i) = 3.*(AL(i)-p(i)) 1603 1656 AR(i) = AL(i) - A6(i) … … 1607 1660 endif 1608 1661 250 continue 1662 END DO 1609 1663 endif 1610 1664 return … … 1617 1671 integer i,j 1618 1672 C 1619 do 35 j=j1,j2 1620 do 35 i=2,IMR 1621 35 CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j)) 1622 C 1623 do 45 j=j1,j2 1624 45 CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j)) 1625 C 1626 do 55 i=1,IMR*JMR 1627 55 CRY(i,2) = DTDY5*(V(i,2)+V(i,1)) 1673 do j=j1,j2 1674 do i=2,IMR 1675 CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j)) 1676 END DO 1677 END DO 1678 C 1679 do j=j1,j2 1680 CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j)) 1681 END DO 1682 C 1683 do i=1,IMR*JMR 1684 CRY(i,2) = DTDY5*(V(i,2)+V(i,1)) 1685 END DO 1628 1686 return 1629 1687 end … … 1636 1694 real ph5 1637 1695 JMR = JNP-1 1638 do 55j=2,JNP1696 do j=2,JNP 1639 1697 ph5 = -0.5*PI + (REAL(J-1)-0.5)*DP 1640 55 cose(j) = cos(ph5) 1698 cose(j) = cos(ph5) 1699 END DO 1641 1700 C 1642 1701 JEQ = (JNP+1) / 2 1643 if(JMR .eq.2*(JMR/2) ) then1702 if(JMR == 2*(JMR/2) ) then 1644 1703 do j=JNP, JEQ+1, -1 1645 1704 cose(j) = cose(JNP+2-j) … … 1653 1712 endif 1654 1713 C 1655 do 66 j=2,JMR 1656 66 cosp(j) = 0.5*(cose(j)+cose(j+1)) 1714 do j=2,JMR 1715 cosp(j) = 0.5*(cose(j)+cose(j+1)) 1716 END DO 1657 1717 cosp(1) = 0. 1658 1718 cosp(JNP) = 0. … … 1668 1728 C 1669 1729 phi = -0.5*PI 1670 do 55j=2,JNP-11730 do j=2,JNP-1 1671 1731 phi = phi + DP 1672 55 cosp(j) = cos(phi) 1732 cosp(j) = cos(phi) 1733 END DO 1673 1734 cosp( 1) = 0. 1674 1735 cosp(JNP) = 0. 1675 1736 C 1676 do 66j=2,JNP1737 do j=2,JNP 1677 1738 cose(j) = 0.5*(cosp(j)+cosp(j-1)) 1678 66 CONTINUE 1679 C 1680 do 77j=2,JNP-11739 END DO 1740 C 1741 do j=2,JNP-1 1681 1742 cosp(j) = 0.5*(cose(j)+cose(j+1)) 1682 77 CONTINUE 1743 END DO 1683 1744 return 1684 1745 end … … 1702 1763 icr = 1 1703 1764 call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) 1704 if(ipy .eq.0) goto 501765 if(ipy==0) goto 50 1705 1766 call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) 1706 if(ipx .eq.0) goto 501767 if(ipx==0) goto 50 1707 1768 C 1708 1769 if(cross) then 1709 1770 call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) 1710 1771 endif 1711 if(icr .eq.0) goto 501772 if(icr==0) goto 50 1712 1773 C 1713 1774 C Vertical filling... 1714 1775 do i=1,len 1715 IF( Q(i,j1,1) .LT.0.) THEN1776 IF( Q(i,j1,1)<0.) THEN 1716 1777 ip = ip + 1 1717 1778 Q(i,j1,2) = Q(i,j1,2) + Q(i,j1,1) … … 1721 1782 C 1722 1783 50 continue 1723 DO 225L = 2,NLAYM11784 DO L = 2,NLAYM1 1724 1785 icr = 1 1725 1786 C 1726 1787 call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) 1727 if(ipy .eq.0) goto 2251788 if(ipy==0) goto 225 1728 1789 call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) 1729 if(ipx .eq.0) go to 2251790 if(ipx==0) go to 225 1730 1791 if(cross) then 1731 1792 call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) 1732 1793 endif 1733 if(icr .eq.0) goto 2251794 if(icr==0) goto 225 1734 1795 C 1735 1796 do i=1,len 1736 IF( Q(I,j1,L) .LT.0.) THEN1797 IF( Q(I,j1,L)<0.) THEN 1737 1798 C 1738 1799 ip = ip + 1 … … 1749 1810 ENDDO 1750 1811 225 CONTINUE 1812 END DO 1751 1813 C 1752 1814 C BOTTOM LAYER … … 1755 1817 C 1756 1818 call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) 1757 if(ipy .eq.0) goto 9111819 if(ipy==0) goto 911 1758 1820 call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) 1759 if(ipx .eq.0) goto 9111821 if(ipx==0) goto 911 1760 1822 C 1761 1823 call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) 1762 if(icr .eq.0) goto 9111824 if(icr==0) goto 911 1763 1825 C 1764 1826 DO I=1,len 1765 IF( Q(I,j1,L) .LT.0.) THEN1827 IF( Q(I,j1,L)<0.) THEN 1766 1828 ip = ip + 1 1767 1829 c … … 1780 1842 911 continue 1781 1843 C 1782 if(ip .gt.IMR) then1844 if(ip>IMR) then 1783 1845 write(6,*) 'IC=',IC,' STEP=',NSTEP, 1784 1846 & ' Vertical filling pts=',ip 1785 1847 endif 1786 1848 C 1787 if(sum .gt.1.e-25) then1849 if(sum>1.e-25) then 1788 1850 write(6,*) IC,NSTEP,' Mass source from the ground=',sum 1789 1851 endif … … 1798 1860 real :: dq,dn,d0,d1,ds,d2 1799 1861 icr = 0 1800 do 65j=j1+1,j2-11801 DO 50i=1,IMR-11802 IF(q(i,j) .LT.0.) THEN1862 do j=j1+1,j2-1 1863 DO i=1,IMR-1 1864 IF(q(i,j)<0.) THEN 1803 1865 icr = 1 1804 1866 dq = - q(i,j)*cosp(j) … … 1816 1878 q(i,j) = (d2 - dq)*acosp(j) + tiny 1817 1879 endif 1818 50 continue 1819 if(icr .eq.0 .and. q(IMR,j).ge.0.) goto 651820 DO 55i=2,IMR1821 IF(q(i,j) .LT.0.) THEN1880 END DO 1881 if(icr==0 .and. q(IMR,j)>=0.) goto 65 1882 DO i=2,IMR 1883 IF(q(i,j)<0.) THEN 1822 1884 icr = 1 1823 1885 dq = - q(i,j)*cosp(j) … … 1835 1897 q(i,j) = (d2 - dq)*acosp(j) + tiny 1836 1898 endif 1837 55 continue 1899 END DO 1838 1900 C ***************************************** 1839 1901 C i=1 1840 1902 i=1 1841 IF(q(i,j) .LT.0.) THEN1903 IF(q(i,j)<0.) THEN 1842 1904 icr = 1 1843 1905 dq = - q(i,j)*cosp(j) … … 1858 1920 C i=IMR 1859 1921 i=IMR 1860 IF(q(i,j) .LT.0.) THEN1922 IF(q(i,j)<0.) THEN 1861 1923 icr = 1 1862 1924 dq = - q(i,j)*cosp(j) … … 1876 1938 C ***************************************** 1877 1939 65 continue 1878 C 1879 do i=1,IMR 1880 if(q(i,j1).lt.0. .or. q(i,j2).lt.0.) then 1940 END DO 1941 C 1942 do i=1,IMR 1943 if(q(i,j1)<0. .or. q(i,j2)<0.) then 1881 1944 icr = 1 1882 1945 goto 80 … … 1886 1949 80 continue 1887 1950 C 1888 if(q(1,1) .lt.0. .or. q(1,jnp).lt.0.) then1951 if(q(1,1)<0. .or. q(1,jnp)<0.) then 1889 1952 icr = 1 1890 1953 endif … … 1910 1973 C 1911 1974 ipy = 0 1912 do 55j=j1+1,j2-11913 DO 55i=1,IMR1914 IF(q(i,j) .LT.0.) THEN1975 do j=j1+1,j2-1 1976 DO i=1,IMR 1977 IF(q(i,j)<0.) THEN 1915 1978 ipy = 1 1916 1979 dq = - q(i,j)*cosp(j) … … 1928 1991 q(i,j) = (d2 - dq)*acosp(j) + tiny 1929 1992 endif 1930 55 continue 1993 END DO 1994 END DO 1931 1995 C 1932 1996 do i=1,imr 1933 IF(q(i,j1) .LT.0.) THEN1997 IF(q(i,j1)<0.) THEN 1934 1998 ipy = 1 1935 1999 dq = - q(i,j1)*cosp(j1) … … 1945 2009 j = j2 1946 2010 do i=1,imr 1947 IF(q(i,j) .LT.0.) THEN2011 IF(q(i,j)<0.) THEN 1948 2012 ipy = 1 1949 2013 dq = - q(i,j)*cosp(j) … … 1958 2022 C 1959 2023 C Check Poles. 1960 if(q(1,1) .lt.0.) then2024 if(q(1,1)<0.) then 1961 2025 dq = q(1,1)*cap1/REAL(IMR)*acosp(j1) 1962 2026 do i=1,imr 1963 2027 q(i,1) = 0. 1964 2028 q(i,j1) = q(i,j1) + dq 1965 if(q(i,j1) .lt.0.) ipy = 11966 enddo 1967 endif 1968 C 1969 if(q(1,JNP) .lt.0.) then2029 if(q(i,j1)<0.) ipy = 1 2030 enddo 2031 endif 2032 C 2033 if(q(1,JNP)<0.) then 1970 2034 dq = q(1,JNP)*cap1/REAL(IMR)*acosp(j2) 1971 2035 do i=1,imr 1972 2036 q(i,JNP) = 0. 1973 2037 q(i,j2) = q(i,j2) + dq 1974 if(q(i,j2) .lt.0.) ipy = 12038 if(q(i,j2)<0.) ipy = 1 1975 2039 enddo 1976 2040 endif … … 1988 2052 ipx = 0 1989 2053 C Copy & swap direction for vectorization. 1990 do 25 i=1,imr 1991 do 25 j=j1,j2 1992 25 qtmp(j,i) = q(i,j) 1993 C 1994 do 55 i=2,imr-1 1995 do 55 j=j1,j2 1996 if(qtmp(j,i).lt.0.) then 2054 do i=1,imr 2055 do j=j1,j2 2056 qtmp(j,i) = q(i,j) 2057 END DO 2058 END DO 2059 C 2060 do i=2,imr-1 2061 do j=j1,j2 2062 if(qtmp(j,i)<0.) then 1997 2063 ipx = 1 1998 2064 c west … … 2007 2073 qtmp(j,i) = qtmp(j,i) + d2 + tiny 2008 2074 endif 2009 55 continue 2075 END DO 2076 END DO 2010 2077 c 2011 2078 i=1 2012 do 65j=j1,j22013 if(qtmp(j,i) .lt.0.) then2079 do j=j1,j2 2080 if(qtmp(j,i)<0.) then 2014 2081 ipx = 1 2015 2082 c west … … 2025 2092 qtmp(j,i) = qtmp(j,i) + d2 + tiny 2026 2093 endif 2027 65 continue 2094 END DO 2028 2095 i=IMR 2029 do 75j=j1,j22030 if(qtmp(j,i) .lt.0.) then2096 do j=j1,j2 2097 if(qtmp(j,i)<0.) then 2031 2098 ipx = 1 2032 2099 c west … … 2042 2109 qtmp(j,i) = qtmp(j,i) + d2 + tiny 2043 2110 endif 2044 75 continue 2045 C 2046 if(ipx.ne.0) then 2047 do 85 j=j1,j2 2048 do 85 i=1,imr 2049 85 q(i,j) = qtmp(j,i) 2111 END DO 2112 C 2113 if(ipx/=0) then 2114 do j=j1,j2 2115 do i=1,imr 2116 q(i,j) = qtmp(j,i) 2117 END DO 2118 END DO 2050 2119 else 2051 2120 C 2052 2121 C Poles. 2053 if(q(1,1) .lt.0. or. q(1,JNP).lt.0.) ipx = 12122 if(q(1,1)<0 .or. q(1,JNP)<0.) ipx = 1 2054 2123 endif 2055 2124 return … … 2065 2134 integer IC,k,i 2066 2135 C 2067 do 4000IC = 1, nc2068 C 2069 do 1000k=1,km2070 do 1000i=1,im2136 do IC = 1, nc 2137 C 2138 do k=1,km 2139 do i=1,im 2071 2140 qtmp(i,k) = q(i,km+1-k,IC) 2072 1000 continue 2073 C 2074 do 2000 i=1,im*km 2075 2000 q(i,1,IC) = qtmp(i,1) 2076 4000 continue 2141 END DO 2142 END DO 2143 C 2144 do i=1,im*km 2145 q(i,1,IC) = qtmp(i,1) 2146 END DO 2147 END DO 2077 2148 return 2078 2149 end -
LMDZ6/trunk/libf/filtrez/inifgn.F
r2598 r5079 28 28 pi = 2.* ASIN(1.) 29 29 C 30 DO 5i=1,iim30 DO i=1,iim 31 31 dlonu(i)= xprimu( i ) 32 32 dlonv(i)= xprimv( i ) 33 5 CONTINUE33 END DO 34 34 35 DO 12i=1,iim35 DO i=1,iim 36 36 sddv(i) = SQRT(dlonv(i)) 37 37 sddu(i) = SQRT(dlonu(i)) 38 38 unsddu(i) = 1./sddu(i) 39 39 unsddv(i) = 1./sddv(i) 40 12 CONTINUE40 END DO 41 41 C 42 DO 17j=1,iim43 DO 17i=1,iim42 DO j=1,iim 43 DO i=1,iim 44 44 vec(i,j) = 0. 45 45 vec1(i,j) = 0. 46 46 eignfnv(i,j) = 0. 47 47 eignfnu(i,j) = 0. 48 17 CONTINUE 48 END DO 49 END DO 49 50 c 50 51 c 51 52 eignfnv(1,1) = -1. 52 53 eignfnv(iim,1) = 1. 53 DO 20i=1,imm154 DO i=1,imm1 54 55 eignfnv(i+1,i+1)= -1. 55 56 eignfnv(i,i+1) = 1. 56 20 CONTINUE57 DO 25j=1,iim58 DO 25i=1,iim57 END DO 58 DO j=1,iim 59 DO i=1,iim 59 60 eignfnv(i,j) = eignfnv(i,j)/(sddu(i)*sddv(j)) 60 25 CONTINUE 61 DO 30 j=1,iim 62 DO 30 i=1,iim 61 END DO 62 END DO 63 DO j=1,iim 64 DO i=1,iim 63 65 eignfnu(i,j) = -eignfnv(j,i) 64 30 CONTINUE 66 END DO 67 END DO 65 68 c 66 69 #ifdef CRAY
Note: See TracChangeset
for help on using the changeset viewer.