Changeset 2270 for LMDZ5/trunk/libf/dyn3dmem
- Timestamp:
- May 7, 2015, 5:45:04 PM (9 years ago)
- Location:
- LMDZ5/trunk/libf/dyn3dmem
- Files:
-
- 1 added
- 12 edited
-
advtrac_loc.F (modified) (4 diffs)
-
caladvtrac_loc.F (modified) (3 diffs)
-
check_isotopes_loc.F (added)
-
dynetat0_loc.F (modified) (2 diffs)
-
iniacademic_loc.F90 (modified) (3 diffs)
-
integrd_loc.F (modified) (8 diffs)
-
leapfrog_loc.F (modified) (22 diffs)
-
qminimum_loc.F (modified) (9 diffs)
-
vlsplt_loc.F (modified) (43 diffs)
-
vlspltgen_loc.F (modified) (12 diffs)
-
vlspltgen_mod.F90 (modified) (2 diffs)
-
vlspltqs_loc.F (modified) (26 diffs)
-
vlz_mod.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/dyn3dmem/advtrac_loc.F
r1987 r2270 24 24 USE Vampir 25 25 USE times 26 USE infotrac, ONLY: nqtot, iadv 26 USE infotrac, ONLY: nqtot, iadv, ok_iso_verif 27 27 USE control_mod, ONLY: iapp_tracvl, day_step, planet_type 28 28 USE advtrac_mod, ONLY: finmasse … … 82 82 !$OMP THREADPRIVATE(testRequest) 83 83 84 c test sur l' eventuelle creation de valeurs negatives de la masse84 c test sur l''eventuelle creation de valeurs negatives de la masse 85 85 ijb=ij_begin 86 86 ije=ij_end … … 155 155 c$OMP BARRIER 156 156 157 write(*,*) 'advtrac 157: appel de vlspltgen_loc' 157 158 call vlspltgen_loc( q,iadv, 2., massem, wg , 158 159 * pbarug,pbarvg,dtvr,p, 159 160 * pk,teta ) 161 162 write(*,*) 'advtrac 162: apres appel vlspltgen_loc' 163 if (ok_iso_verif) then 164 call check_isotopes(q,ijb_u,ije_u,'advtrac 162') 165 endif !if (ok_iso_verif) then 160 166 161 167 #ifdef DEBUG_IO … … 356 362 c$OMP END DO 357 363 358 CALL qminimum_loc( q, 2, finmasse ) 364 ! CRisi: on passe nqtot et non nq 365 CALL qminimum_loc( q, nqtot, finmasse ) 359 366 360 367 endif ! of if (planet_type=="earth") -
LMDZ5/trunk/libf/dyn3dmem/caladvtrac_loc.F
r1907 r2270 56 56 !$OMP THREADPRIVATE(Request_vanleer) 57 57 58 58 write(*,*) 'caladvtrac 58: entree' 59 59 ijbu=ij_begin 60 60 ijeu=ij_end … … 109 109 110 110 IF ( iadvtr.EQ.iapp_tracvl ) THEN 111 write(*,*) 'caladvtrac 133' 111 112 c$OMP MASTER 112 113 call suspend_timer(timer_caldyn) … … 183 184 CALL WriteField_u('wg1',wg_adv) 184 185 #endif 186 write(*,*) 'caladvtrac 185' 185 187 CALL advtrac_loc( pbarug_adv,pbarvg_adv,wg_adv, 186 188 * p_adv, massem_adv,q_adv, teta_adv, 187 . pk_adv) 189 . pk_adv) 190 write(*,*) 'caladvtrac 189' 188 191 189 192 -
LMDZ5/trunk/libf/dyn3dmem/dynetat0_loc.F
r1939 r2270 366 366 write(lunout,*)"Il est donc initialise a zero" 367 367 q(:,:,iq)=0. 368 369 ! CRisi: pour les isotopes, on peut faire init théorique 370 ! distill de Rayleigh très simplifiée 371 if (ok_isotopes) then 372 if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.0)) then 373 q(:,:,iq)=q(:,:,iqpere(iq)) & 374 & *tnat(iso_num(iq)) & 375 & *(q(:,:,iqpere(iq))/30.e-3) & 376 & **(alpha_ideal(iso_num(iq))-1) 377 endif 378 if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.1)) then 379 q(:,:,iq)=q(:,:,iqiso(iso_indnum(iq),phase_num(iq))) 380 endif 381 endif !if (ok_isotopes) then 382 368 383 ELSE 369 384 #ifdef NC_DOUBLE … … 380 395 381 396 ENDIF 382 ENDDO 397 ENDDO !DO iq=1,nqtot 398 399 if (ok_iso_verif) then 400 call check_isotopes(q,ijb_u,ije_u,'dynetat0_loc') 401 endif !if (ok_iso_verif) then 383 402 384 403 DEALLOCATE(q_glo) -
LMDZ5/trunk/libf/dyn3dmem/iniacademic_loc.F90
r2087 r2270 7 7 use exner_hyb_m, only: exner_hyb 8 8 use exner_milieu_m, only: exner_milieu 9 USE infotrac, ONLY : nqtot 9 USE infotrac, ONLY: nqtot,niso_possibles,ok_isotopes,iqpere,ok_iso_verif,tnat,alpha_ideal, & 10 & iqiso,phase_num,iso_indnum,iso_num,zone_num 10 11 USE control_mod, ONLY: day_step,planet_type 11 12 USE parallel_lmdz, ONLY: ijb_u, ije_u, ijb_v, ije_v … … 110 111 ztot0 = 0. 111 112 stot0 = 0. 112 ang0 = 0. 113 ang0 = 0. 113 114 114 115 if (llm == 1) then … … 269 270 if (planet_type=="earth") then 270 271 ! Earth: first two tracers will be water 272 271 273 do i=1,nqtot 272 274 if (i == 1) q(ijb_u:ije_u,:,i)=1.e-10 273 275 if (i == 2) q(ijb_u:ije_u,:,i)=1.e-15 274 276 if (i.gt.2) q(ijb_u:ije_u,:,i)=0. 277 278 ! CRisi: init des isotopes 279 ! distill de Rayleigh très simplifiée 280 if (ok_isotopes) then 281 if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then 282 q(ijb_u:ije_u,:,i)=q(ijb_u:ije_u,:,iqpere(i)) & 283 & *tnat(iso_num(i)) & 284 & *(q(ijb_u:ije_u,:,iqpere(i))/30.e-3) & 285 & **(alpha_ideal(iso_num(i))-1) 286 endif 287 if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then 288 q(ijb_u:ije_u,:,i)=q(ijb_u:ije_u,:,iqiso(iso_indnum(i),phase_num(i))) 289 endif 290 endif !if (ok_isotopes) then 291 275 292 enddo 276 293 else 277 294 q(ijb_u:ije_u,:,:)=0 278 295 endif ! of if (planet_type=="earth") 296 297 if (ok_iso_verif) then 298 call check_isotopes(q,ijb_u,ije_u,'iniacademic_loc') 299 endif !if (ok_iso_verif) then 279 300 280 301 ! add random perturbation to temperature -
LMDZ5/trunk/libf/dyn3dmem/integrd_loc.F
r2110 r2270 11 11 USE write_field 12 12 USE integrd_mod 13 USE infotrac, ONLY: ok_iso_verif ! ajout CRisi 13 14 IMPLICIT NONE 14 15 … … 86 87 INTEGER :: ierr 87 88 89 !write(*,*) 'integrd 88: entree, nq=',nq 88 90 c----------------------------------------------------------------------- 91 92 if (ok_iso_verif) then 93 call check_isotopes(q,ijb,ije,'integrd 92') 94 endif !if (ok_iso_verif) then 95 89 96 c$OMP BARRIER 90 97 if (pole_nord) THEN … … 125 132 DO 2 ij = ijb,ije 126 133 pscr (ij) = ps0(ij) 127 ps (ij) = psm1(ij) + dt * dp(ij) 134 ps (ij) = psm1(ij) + dt * dp(ij) 135 128 136 2 CONTINUE 137 129 138 c$OMP END DO 130 139 c$OMP BARRIER … … 214 223 c$OMP END MASTER 215 224 c$OMP BARRIER 225 !write(*,*) 'integrd 217' 216 226 c 217 227 c ... Calcul de la nouvelle masse d'air au dernier temps integre t+1 ... … … 219 229 220 230 CALL pression_loc ( ip1jmp1, ap, bp, ps, p ) 231 221 232 c$OMP BARRIER 222 233 CALL massdair_loc ( p , masse ) … … 334 345 ENDDO 335 346 ENDDO 347 336 348 c$OMP END DO NOWAIT 337 349 c$OMP BARRIER 338 350 351 if (ok_iso_verif) then 352 call check_isotopes(q,ijb,ije,'integrd 342') 353 endif !if (ok_iso_verif) then 354 355 !write(*,*) 'integrd 341' 339 356 CALL qminimum_loc( q, nq, deltap ) 357 !write(*,*) 'integrd 343' 358 359 if (ok_iso_verif) then 360 call check_isotopes(q,ijb,ije,'integrd 346') 361 endif !if (ok_iso_verif) then 340 362 c 341 363 c ..... Calcul de la valeur moyenne, unique aux poles pour q ..... … … 387 409 388 410 ENDIF 411 412 if (ok_iso_verif) then 413 call check_isotopes(q,ijb,ije,'integrd 409') 414 endif !if (ok_iso_verif) then 389 415 390 416 ! Ehouarn: forget about finvmaold … … 404 430 405 431 15 continue 432 !write(*,*) 'integrd 410' 406 433 407 434 c$OMP DO SCHEDULE(STATIC) -
LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F
r2221 r2270 200 200 LOGICAL, SAVE :: firstcall=.TRUE. 201 201 TYPE(distrib),SAVE :: new_dist 202 203 if (ok_iso_verif) then 204 call check_isotopes(q0,ijb_u,ije_u,'leapfrog204: debut') 205 endif !if (ok_iso_verif) then 202 206 203 207 c$OMP MASTER … … 219 223 itaufinp1 = itaufin +1 220 224 225 if (ok_iso_verif) then 226 call check_isotopes(q0,ijb_u,ije_u,'leapfrog 226') 227 endif !if (ok_iso_verif) then 228 221 229 itau = 0 222 230 physic=.true. … … 231 239 phis=phis0 232 240 q=q0 241 242 if (ok_iso_verif) then 243 call check_isotopes(q,ijb_u,ije_u,'leapfrog 239') 244 endif !if (ok_iso_verif) then 233 245 234 246 ! iday = day_ini+itau/day_step … … 296 308 297 309 1 CONTINUE ! Matsuno Forward step begins here 298 310 write(*,*) 'leapfrog 298: itau=',itau 299 311 jD_cur = jD_ref + day_ini - day_ref + & 300 312 & itau/day_step … … 306 318 endif 307 319 320 if (ok_iso_verif) then 321 call check_isotopes(q,ijb_u,ije_u,'leapfrog 321') 322 endif !if (ok_iso_verif) then 308 323 309 324 #ifdef CPP_IOIPSL … … 384 399 cym call minmax(ijp1llm,q(:,:,3),zqmin,zqmax) 385 400 401 402 if (ok_iso_verif) then 403 call check_isotopes(q,ijb_u,ije_u,'leapfrog 400') 404 endif !if (ok_iso_verif) then 405 386 406 2 CONTINUE ! Matsuno backward or leapfrog step begins here 407 408 409 if (ok_iso_verif) then 410 call check_isotopes(q,ijb_u,ije_u,'leapfrog 402') 411 endif !if (ok_iso_verif) then 387 412 388 413 c$OMP MASTER … … 455 480 c$OMP END MASTER 456 481 482 483 if (ok_iso_verif) then 484 call check_isotopes(q,ijb_u,ije_u,'leapfrog 471') 485 endif !if (ok_iso_verif) then 457 486 458 487 !ym PAS D'AJUSTEMENT POUR LE MOMENT … … 574 603 575 604 605 if (ok_iso_verif) then 606 call check_isotopes(q,ijb_u,ije_u,'leapfrog 589') 607 endif !if (ok_iso_verif) then 576 608 577 609 c----------------------------------------------------------------------- … … 635 667 ! compute geopotential phi() 636 668 CALL geopot_loc ( ip1jmp1, teta , pk , pks, phis , phi ) 637 669 670 if (ok_iso_verif) then 671 call check_isotopes(q,ijb_u,ije_u,'leapfrog 651') 672 endif !if (ok_iso_verif) then 638 673 639 674 call VTb(VTcaldyn) … … 644 679 ! CALL FTRACE_REGION_BEGIN("caldyn") 645 680 time = jD_cur + jH_cur 681 646 682 CALL caldyn_loc 647 683 $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis , … … 670 706 c ------------------------------------------------------------- 671 707 708 if (ok_iso_verif) then 709 call check_isotopes(q,ijb_u,ije_u, 710 & 'leapfrog 686: avant caladvtrac') 711 endif !if (ok_iso_verif) then 672 712 673 713 IF( forward. OR . leapf ) THEN 674 714 ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step 715 write(*,*) 'leapfrog 679: avant CALL caladvtrac_loc' 675 716 CALL caladvtrac_loc(q,pbaru,pbarv, 676 717 * p, masse, dq, teta, 677 718 . flxw,pk, iapptrac) 719 720 write(*,*) 'leapfrog 719' 721 if (ok_iso_verif) then 722 call check_isotopes(q,ijb_u,ije_u, 723 & 'leapfrog 698: apres caladvtrac') 724 endif !if (ok_iso_verif) then 678 725 679 726 ! do j=1,nqtot … … 708 755 ! CALL FTRACE_REGION_BEGIN("integrd") 709 756 710 CALL integrd_loc ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 , 757 write(*,*) 'leapfrog 720' 758 if (ok_iso_verif) then 759 call check_isotopes(q,ijb_u,ije_u,'leapfrog 756') 760 endif !if (ok_iso_verif) then 761 762 ! CRisi: pourquoi aller jusqu'à 2 et non pas jusqu'à nqtot?? 763 CALL integrd_loc ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 , 711 764 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis) 712 765 ! $ finvmaold ) 713 766 767 write(*,*) 'leapfrog 724' 768 if (ok_iso_verif) then 769 call check_isotopes(q,ijb_u,ije_u,'leapfrog 762') 770 endif !if (ok_iso_verif) then 771 714 772 ! CALL FTRACE_REGION_END("integrd") 715 773 c$OMP BARRIER … … 724 782 call WriteField_u('ps_int',ps) 725 783 #endif 784 785 if (ok_iso_verif) then 786 call check_isotopes(q,ijb_u,ije_u,'leapfrog 775') 787 endif !if (ok_iso_verif) then 788 726 789 c do j=1,nqtot 727 790 c call WriteField_p('q'//trim(int2str(j)), … … 1082 1145 ENDIF ! of IF( apphys ) 1083 1146 1147 if (ok_iso_verif) then 1148 call check_isotopes(q,ijb_u,ije_u,'leapfrog 1132') 1149 endif !if (ok_iso_verif) then 1150 write(*,*) 'leapfrog 1134: iflag_phys=',iflag_phys 1151 1084 1152 IF(iflag_phys.EQ.2) THEN ! "Newtonian" case 1085 1153 c$OMP MASTER … … 1146 1214 1147 1215 cc$OMP END PARALLEL 1216 if (ok_iso_verif) then 1217 call check_isotopes(q,ijb_u,ije_u,'leapfrog 1196') 1218 endif !if (ok_iso_verif) then 1148 1219 1149 1220 c----------------------------------------------------------------------- 1150 1221 c dissipation horizontale et verticale des petites echelles: 1151 1222 c ---------------------------------------------------------- 1152 1223 write(*,*) 'leapfrog 1163: apdiss=',apdiss 1153 1224 IF(apdiss) THEN 1154 1225 … … 1379 1450 c call abort_gcm(modname,abort_message,0) 1380 1451 c ENDIF 1381 1452 1453 if (ok_iso_verif) then 1454 call check_isotopes(q,ijb_u,ije_u,'leapfrog 1430') 1455 endif !if (ok_iso_verif) then 1456 1382 1457 c ******************************************************************** 1383 1458 c ******************************************************************** … … 1455 1530 ENDIF 1456 1531 1532 if (ok_iso_verif) then 1533 call check_isotopes(q,ijb_u,ije_u,'leapfrog 1509') 1534 endif !if (ok_iso_verif) then 1535 1457 1536 IF ( .NOT.purmats ) THEN 1458 1537 c ........................................................ … … 1526 1605 ENDIF 1527 1606 1607 if (ok_iso_verif) then 1608 call check_isotopes(q,ijb_u,ije_u,'leapfrog 1584') 1609 endif !if (ok_iso_verif) then 1610 1528 1611 c----------------------------------------------------------------------- 1529 1612 c ecriture de la bande histoire: … … 1562 1645 ENDIF ! of IF (itau.EQ.itaufin) 1563 1646 1647 if (ok_iso_verif) then 1648 call check_isotopes(q,ijb_u,ije_u,'leapfrog 1624') 1649 endif !if (ok_iso_verif) then 1650 1564 1651 c----------------------------------------------------------------------- 1565 1652 c gestion de l'integration temporelle: … … 1596 1683 1597 1684 ELSE ! of IF (.not.purmats) 1685 1686 1687 if (ok_iso_verif) then 1688 call check_isotopes(q,ijb_u,ije_u,'leapfrog 1664') 1689 endif !if (ok_iso_verif) then 1598 1690 1599 1691 c ........................................................ … … 1631 1723 1632 1724 ELSE ! of IF(forward) i.e. backward step 1725 1726 1727 if (ok_iso_verif) then 1728 call check_isotopes(q,ijb_u,ije_u,'leapfrog 1698') 1729 endif !if (ok_iso_verif) then 1633 1730 1634 1731 IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN … … 1683 1780 ENDIF ! of IF (forward) 1684 1781 1782 1783 if (ok_iso_verif) then 1784 call check_isotopes(q,ijb_u,ije_u,'leapfrog 1750') 1785 endif !if (ok_iso_verif) then 1786 1685 1787 END IF ! of IF(.not.purmats) 1686 1788 c$OMP MASTER -
LMDZ5/trunk/libf/dyn3dmem/qminimum_loc.F
r1907 r2270 1 SUBROUTINE qminimum_loc( q,nq ,deltap )1 SUBROUTINE qminimum_loc( q,nqtot,deltap ) 2 2 USE parallel_lmdz 3 USE infotrac, ONLY: ok_isotopes,ntraciso,iqiso,ok_iso_verif 3 4 IMPLICIT none 4 5 c … … 10 11 #include "comvert.h" 11 12 c 12 INTEGER nq 13 REAL q(ijb_u:ije_u,llm,nq ), deltap(ijb_u:ije_u,llm)13 INTEGER nqtot ! CRisi: on remplace nq par nqtot 14 REAL q(ijb_u:ije_u,llm,nqtot), deltap(ijb_u:ije_u,llm) 14 15 c 15 16 INTEGER iq_vap, iq_liq … … 27 28 INTEGER i, k, iq 28 29 REAL zx_defau, zx_abc, zx_pump(ijb_u:ije_u), pompe 30 31 real zx_defau_diag(ijb_u:ije_u,llm,2) 32 real q_follow(ijb_u:ije_u,llm,2) 29 33 c 30 34 REAL SSUM … … 38 42 INTEGER Index_pump(ij_end-ij_begin+1) 39 43 INTEGER nb_pump 44 INTEGER ixt 45 INTEGER iso_verif_noNaN_nostop 40 46 c 41 47 c Quand l'eau liquide est trop petite (ou negative), on prend … … 44 50 c 45 51 52 write(*,*) 'qminimum 52: entree' 53 if (ok_iso_verif) then 54 call check_isotopes(q,ij_begin,ij_end,'qminimum 52') 55 endif !if (ok_iso_verif) then 56 46 57 ijb=ij_begin 47 58 ije=ij_end 48 59 60 zx_defau_diag(ijb:ije,:,:)=0.0 61 q_follow(ijb:ije,:,:)=q(ijb:ije,:,:) 62 63 !write(*,*) 'qminimum 57' 49 64 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 50 65 DO 1000 k = 1, llm 51 66 DO 1040 i = ijb, ije 52 67 if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then 68 69 if (ok_isotopes) then 70 zx_defau_diag(i,k,iq_liq)=AMAX1 71 : ( seuil_liq - q(i,k,iq_liq), 0.0 ) 72 endif !if (ok_isotopes) then 73 53 74 q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq 54 75 q(i,k,iq_liq) = seuil_liq … … 60 81 c ---> SYNCHRO OPENMP ICI 61 82 83 62 84 c 63 85 c Quand l'eau vapeur est trop faible (ou negative), on complete 64 86 c le defaut en prennant de l'eau vapeur de la couche au-dessous. 65 87 c 88 !write(*,*) 'qminimum 81' 66 89 iq = iq_vap 67 90 c … … 70 93 c$OMP DO SCHEDULE(STATIC) 71 94 DO i = ijb, ije 95 72 96 if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then 97 98 if (ok_isotopes) then 99 zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 0.0 ) 100 endif !if (ok_isotopes) then 101 73 102 q(i,k-1,iq) = q(i,k-1,iq) - ( seuil_vap - q(i,k,iq) ) * 74 103 & deltap(i,k) / deltap(i,k-1) 75 104 q(i,k,iq) = seuil_vap 105 76 106 endif 77 107 ENDDO … … 79 109 ENDDO 80 110 c$OMP BARRIER 111 81 112 c 82 113 c Quand il s'agit de la premiere couche au-dessus du sol, on 83 114 c doit imprimer un message d'avertissement (saturation possible). 84 115 c 116 !write(*,*) 'qminimum 106' 85 117 nb_pump=0 86 118 c$OMP DO SCHEDULE(STATIC) … … 103 135 ENDDO 104 136 ENDIF 137 138 !write(*,*) 'qminimum 128' 139 if (ok_isotopes) then 140 ! CRisi: traiter de même les traceurs d'eau 141 ! Mais il faut les prendre à l'envers pour essayer de conserver la 142 ! masse. 143 ! 1) pompage dans le sol 144 ! On suppose que ce pompage se fait sans isotopes -> on ne modifie 145 ! rien ici et on croise les doigts pour que ça ne soit pas trop 146 ! génant 147 DO i = ijb, ije 148 if (zx_pump(i).gt.0.0) then 149 q_follow(i,1,iq_vap)=q_follow(i,1,iq_vap)+zx_pump(i) 150 endif !if (zx_pump(i).gt.0.0) then 151 enddo !DO i = ijb, ije 152 153 ! 2) transfert de vap vers les couches plus hautes 154 !write(*,*) 'qminimum 139' 155 do k=2,llm 156 DO i = ijb, ije 157 if (zx_defau_diag(i,k,iq_vap).gt.0.0) then 158 ! on ajoute la vapeur en k 159 do ixt=1,ntraciso 160 q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap)) 161 : +zx_defau_diag(i,k,iq_vap) 162 : *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap) 163 164 if (ok_iso_verif) then 165 if (iso_verif_noNaN_nostop(q(i,k,iqiso(ixt,iq_vap)), 166 : 'qminimum 155').eq.1) then 167 write(*,*) 'i,k,ixt=',i,k,ixt 168 write(*,*) 'q_follow(i,k-1,iq_vap)=', 169 : q_follow(i,k-1,iq_vap) 170 write(*,*) 'q(i,k,iqiso(ixt,iq_vap))=', 171 : q(i,k,iqiso(ixt,iq_vap)) 172 write(*,*) 'zx_defau_diag(i,k,iq_vap)=', 173 : zx_defau_diag(i,k,iq_vap) 174 write(*,*) 'q(i,k-1,iqiso(ixt,iq_vap))=', 175 : q(i,k-1,iqiso(ixt,iq_vap)) 176 stop 177 endif 178 endif 179 180 ! et on la retranche en k-1 181 q(i,k-1,iqiso(ixt,iq_vap))=q(i,k-1,iqiso(ixt,iq_vap)) 182 : -zx_defau_diag(i,k,iq_vap) 183 : *deltap(i,k)/deltap(i,k-1) 184 : *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap) 185 186 if (ok_iso_verif) then 187 if (iso_verif_noNaN_nostop(q(i,k-1,iqiso(ixt,iq_vap)), 188 : 'qminimum 175').eq.1) then 189 write(*,*) 'k,i,ixt=',k,i,ixt 190 write(*,*) 'q_follow(i,k-1,iq_vap)=', 191 : q_follow(i,k-1,iq_vap) 192 write(*,*) 'q(i,k,iqiso(ixt,iq_vap))=', 193 : q(i,k,iqiso(ixt,iq_vap)) 194 write(*,*) 'zx_defau_diag(i,k,iq_vap)=', 195 : zx_defau_diag(i,k,iq_vap) 196 write(*,*) 'q(i,k-1,iqiso(ixt,iq_vap))=', 197 : q(i,k-1,iqiso(ixt,iq_vap)) 198 stop 199 endif 200 endif 201 202 enddo !do ixt=1,niso 203 q_follow(i,k,iq_vap)= q_follow(i,k,iq_vap) 204 : +zx_defau_diag(i,k,iq_vap) 205 q_follow(i,k-1,iq_vap)= q_follow(i,k-1,iq_vap) 206 : -zx_defau_diag(i,k,iq_vap) 207 : *deltap(i,k)/deltap(i,k-1) 208 endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then 209 enddo !DO i = 1, ip1jmp1 210 enddo !do k=2,llm 211 212 if (ok_iso_verif) then 213 call check_isotopes(q,ijb,ije,'qminimum 168') 214 endif !if (ok_iso_verif) then 215 216 217 ! 3) transfert d'eau de la vapeur au liquide 218 !write(*,*) 'qminimum 164' 219 do k=1,llm 220 DO i = ijb, ije 221 if (zx_defau_diag(i,k,iq_liq).gt.0.0) then 222 223 ! on ajoute eau liquide en k en k 224 do ixt=1,ntraciso 225 q(i,k,iqiso(ixt,iq_liq))=q(i,k,iqiso(ixt,iq_liq)) 226 : +zx_defau_diag(i,k,iq_liq) 227 : *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap) 228 ! et on la retranche à la vapeur en k 229 q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap)) 230 : -zx_defau_diag(i,k,iq_liq) 231 : *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap) 232 enddo !do ixt=1,niso 233 q_follow(i,k,iq_liq)= q_follow(i,k,iq_liq) 234 : +zx_defau_diag(i,k,iq_liq) 235 q_follow(i,k,iq_vap)= q_follow(i,k,iq_vap) 236 : -zx_defau_diag(i,k,iq_liq) 237 endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then 238 enddo !DO i = 1, ip1jmp1 239 enddo !do k=2,llm 240 241 if (ok_iso_verif) then 242 call check_isotopes(q,ijb,ije,'qminimum 197') 243 endif !if (ok_iso_verif) then 244 245 endif !if (ok_isotopes) then 246 !write(*,*) 'qminimum 188' 105 247 c 106 248 RETURN -
LMDZ5/trunk/libf/dyn3dmem/vlsplt_loc.F
r1907 r2270 2 2 ! $Id$ 3 3 ! 4 SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x )4 SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x,iq) 5 5 6 6 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 14 14 c -------------------------------------------------------------------- 15 15 USE parallel_lmdz 16 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 16 17 IMPLICIT NONE 17 18 c … … 25 26 c Arguments: 26 27 c ---------- 27 REAL masse(ijb_u:ije_u,llm),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) 30 REAL w(ijb_u:ije_u,llm) 28 REAL masse(ijb_u:ije_u,llm,nqtot),pente_max 29 REAL u_m( ijb_u:ije_u,llm),pbarv( iip1,jjb_v:jje_v,llm) 30 REAL q(ijb_u:ije_u,llm,nqtot) ! CRisi: ajout dimension nqtot 31 REAL w(ijb_u:ije_u,llm) 32 INTEGER iq ! CRisi 31 33 c 32 34 c Local … … 42 44 REAL u_mq(ijb_u:ije_u,llm) 43 45 46 REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi 47 INTEGER ifils,iq2 ! CRisi 48 44 49 Logical extremum 45 50 … … 51 56 INTEGER ijb,ije,ijb_x,ije_x 52 57 58 write(*,*) 'vlsplt 58: entree dans vlx_loc, iq,ijb_x=', 59 & iq,ijb_x 53 60 c calcul de la pente a droite et a gauche de la maille 54 61 … … 64 71 c calcul des pentes avec limitation, Van Leer scheme I: 65 72 c ----------------------------------------------------- 66 73 ! on a besoin de q entre ijb et ije 67 74 c calcul de la pente aux points u 68 75 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) … … 70 77 71 78 DO ij=ijb,ije-1 72 dxqu(ij)=q(ij+1,l )-q(ij,l)79 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 73 80 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0' 74 c sigu(ij)=u_m(ij,l)/masse(ij,l )81 c sigu(ij)=u_m(ij,l)/masse(ij,l,iq) 75 82 ENDDO 76 83 DO ij=ijb+iip1-1,ije,iip1 … … 126 133 DO l = 1, llm 127 134 DO ij=ijb,ije-1 128 dxqu(ij)=q(ij+1,l )-q(ij,l)135 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 129 136 ENDDO 130 137 DO ij=ijb+iip1-1,ije,iip1 … … 147 154 ENDIF ! (pente_max.lt.-1.e-5) 148 155 156 write(*,*) 'vlx 156: iq,ijb_x=',iq,ijb_x 157 149 158 c bouclage de la pente en iip1: 150 159 c ----------------------------- … … 168 177 DO l=1,llm 169 178 DO ij=ijb,ije-1 170 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l ),171 , 1.+u_m(ij,l)/masse(ij+1,l ),172 , u_m(ij,l ))179 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq), 180 , 1.+u_m(ij,l)/masse(ij+1,l,iq), 181 , u_m(ij,l,iq)) 173 182 zdum(ij,l)=0.5*zdum(ij,l) 174 183 u_mq(ij,l)=cvmgp( 175 , q(ij,l )+zdum(ij,l)*dxq(ij,l),176 , q(ij+1,l )-zdum(ij,l)*dxq(ij+1,l),184 , q(ij,l,iq)+zdum(ij,l)*dxq(ij,l), 185 , q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l), 177 186 , u_m(ij,l)) 178 187 u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l) … … 185 194 c print*,'Cumule ....' 186 195 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 196 ! on a besoin de masse entre ijb et ije 187 197 DO l=1,llm 188 198 DO ij=ijb,ije-1 189 c print*,'masse(',ij,')=',masse(ij,l )199 c print*,'masse(',ij,')=',masse(ij,l,iq) 190 200 IF (u_m(ij,l).gt.0.) THEN 191 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l) 192 u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l)) 201 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq) 202 u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq) 203 : +0.5*zdum(ij,l)*dxq(ij,l)) 193 204 ELSE 194 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l) 195 u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l)) 205 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq) 206 u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq) 207 : -0.5*zdum(ij,l)*dxq(ij+1,l)) 196 208 ENDIF 197 209 ENDDO … … 215 227 c$OMP END DO NOWAIT 216 228 c print*,'Ok test 1' 229 217 230 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 218 231 DO l=1,llm … … 223 236 c$OMP END DO NOWAIT 224 237 c print*,'Ok test 2' 225 238 226 239 227 240 c traitement special pour le cas ou on advecte en longitude plus que le … … 247 260 c & ,'contenu de la maille : ',n0 248 261 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 262 263 249 264 DO l=1,llm 250 265 IF(nl(l).gt.0) THEN … … 258 273 ENDDO 259 274 niju=iju 260 c PRINT*,'niju,nl',niju,nl(l)275 PRINT*,'vlx 278, niju,nl',niju,nl(l) 261 276 262 277 c traitement des mailles … … 270 285 i=ijq-(j-1)*iip1 271 286 c accumulation pour les mailles completements advectees 272 do while(zu_m.gt.masse(ijq,l)) 273 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l) 274 zu_m=zu_m-masse(ijq,l) 287 do while(zu_m.gt.masse(ijq,l,iq)) 288 u_mq(ij,l)=u_mq(ij,l) 289 & +q(ijq,l,iq)*masse(ijq,l,iq) 290 zu_m=zu_m-masse(ijq,l,iq) 275 291 i=mod(i-2+iim,iim)+1 276 292 ijq=(j-1)*iip1+i … … 278 294 c ajout de la maille non completement advectee 279 295 u_mq(ij,l)=u_mq(ij,l)+zu_m* 280 & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l)) 296 & (q(ijq,l,iq)+0.5* 297 & (1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 281 298 ELSE 282 299 ijq=ij+1 283 300 i=ijq-(j-1)*iip1 284 301 c accumulation pour les mailles completements advectees 285 do while(-zu_m.gt.masse(ijq,l)) 286 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l) 287 zu_m=zu_m+masse(ijq,l) 302 do while(-zu_m.gt.masse(ijq,l,iq)) 303 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq) 304 & *masse(ijq,l,iq) 305 zu_m=zu_m+masse(ijq,l,iq) 288 306 i=mod(i,iim)+1 289 307 ijq=(j-1)*iip1+i 290 308 ENDDO 291 309 c ajout de la maille non completement advectee 292 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l )-293 & 0.5*(1.+zu_m/masse(ijq,l ))*dxq(ijq,l))310 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)- 311 & 0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 294 312 ENDIF 295 313 ENDDO … … 299 317 cym ENDIF ! n0.gt.0 300 318 9999 continue 301 302 319 303 320 c bouclage en latitude … … 311 328 c$OMP END DO NOWAIT 312 329 330 ! CRisi: appel récursif de l'advection sur les fils. 331 ! Il faut faire ça avant d'avoir mis à jour q et masse 332 333 write(*,*) 'vlsplt 326: iq,ijb_x,nqfils(iq)=',iq,ijb_x,nqfils(iq) 334 335 if (nqfils(iq).gt.0) then 336 do ifils=1,nqdesc(iq) 337 iq2=iqfils(ifils,iq) 338 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 339 DO l=1,llm 340 DO ij=ijb,ije 341 ! On a besoin de q et masse seulement entre ijb et ije. On ne 342 ! les calcule donc que de ijb à ije 343 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 344 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 345 enddo 346 enddo 347 c$OMP END DO NOWAIT 348 enddo !do ifils=1,nqdesc(iq) 349 do ifils=1,nqfils(iq) 350 iq2=iqfils(ifils,iq) 351 call vlx_loc(Ratio,pente_max,masseq,u_mq,ijb_x,ije_x,iq2) 352 enddo !do ifils=1,nqfils(iq) 353 endif !if (nqfils(iq).gt.0) then 354 ! end CRisi 355 356 write(*,*) 'vlsplt 360: iq,ijb_x=',iq,ijb_x 357 313 358 c calcul des tENDances 314 359 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 315 360 DO l=1,llm 316 361 DO ij=ijb+1,ije 317 new_m=masse(ij,l )+u_m(ij-1,l)-u_m(ij,l)318 q(ij,l )=(q(ij,l)*masse(ij,l)+319 & u_mq(ij-1,l)-u_mq(ij,l))320 & /new_m321 masse(ij,l )=new_m362 new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l) 363 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+ 364 & u_mq(ij-1,l)-u_mq(ij,l)) 365 & /new_m 366 masse(ij,l,iq)=new_m 322 367 ENDDO 323 368 c ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous) 324 369 DO ij=ijb+iip1-1,ije,iip1 325 q(ij-iim,l)=q(ij,l) 326 masse(ij-iim,l)=masse(ij,l) 327 ENDDO 328 ENDDO 329 c$OMP END DO NOWAIT 370 q(ij-iim,l,iq)=q(ij,l,iq) 371 masse(ij-iim,l,iq)=masse(ij,l,iq) 372 ENDDO 373 ENDDO 374 c$OMP END DO NOWAIT 375 write(*,*) 'vlsplt 380: iq,ijb_x=',iq,ijb_x 376 377 ! retablir les fils en rapport de melange par rapport a l'air: 378 ! On calcule q entre ijb+1 et ije -> on fait pareil pour ratio 379 ! puis on boucle en longitude 380 if (nqfils(iq).gt.0) then 381 do ifils=1,nqdesc(iq) 382 iq2=iqfils(ifils,iq) 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 ! DO ij=ijb+iip1-1,ije,iip1 391 enddo !DO l=1,llm 392 c$OMP END DO NOWAIT 393 enddo !do ifils=1,nqdesc(iq) 394 endif !if (nqfils(iq).gt.0) then 395 396 write(*,*) 'vlsplt 399: iq,ijb_x=',iq,ijb_x 330 397 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) 331 398 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1) … … 336 403 337 404 338 SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v )405 SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v,iq) 339 406 c 340 407 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 349 416 c -------------------------------------------------------------------- 350 417 USE parallel_lmdz 418 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 351 419 IMPLICIT NONE 352 420 c … … 361 429 c Arguments: 362 430 c ---------- 363 REAL masse(ijb_u:ije_u,llm ),pente_max431 REAL masse(ijb_u:ije_u,llm,nqtot),pente_max 364 432 REAL masse_adv_v( ijb_v:ije_v,llm) 365 REAL q(ijb_u:ije_u,llm), dq( ijb_u:ije_u,llm) 433 REAL q(ijb_u:ije_u,llm,nqtot), dq( ijb_u:ije_u,llm) 434 INTEGER iq ! CRisi 366 435 c 367 436 c Local … … 392 461 SAVE airej2,airejjm 393 462 c$OMP THREADPRIVATE(airej2,airejjm) 463 464 REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi 465 INTEGER ifils,iq2 ! CRisi 394 466 c 395 467 c … … 401 473 INTEGER ijb,ije 402 474 475 ijb=ij_begin-2*iip1 476 ije=ij_end+2*iip1 477 if (pole_nord) ijb=ij_begin 478 if (pole_sud) ije=ij_end 479 403 480 IF(first) THEN 404 cPRINT*,'Shema Amont nouveau appele dans Vanleer '481 PRINT*,'Shema Amont nouveau appele dans Vanleer ' 405 482 first=.false. 406 483 do i=2,iip1 … … 434 511 if (pole_nord) then 435 512 DO i = 1, iim 436 airescb(i) = aire(i+ iip1) * q(i+ iip1,l )513 airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq) 437 514 ENDDO 438 515 qpns = SSUM( iim, airescb ,1 ) / airej2 … … 441 518 if (pole_sud) then 442 519 DO i = 1, iim 443 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l )520 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq) 444 521 ENDDO 445 522 qpsn = SSUM( iim, airesch ,1 ) / airejjm 446 523 endif 447 524 448 449 450 525 c calcul des pentes aux points v 451 526 … … 455 530 if (pole_sud) ije=ij_end-iip1 456 531 532 ! on a besoin de q entre ij_begin-2*iip1 et ij_end+2*iip1 533 ! Si pole sud, entre ij_begin-2*iip1 et ij_end 534 ! Si pole Nord, entre ij_begin et ij_end+2*iip1 457 535 DO ij=ijb,ije 458 dyqv(ij)=q(ij,l )-q(ij+iip1,l)536 dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq) 459 537 adyqv(ij)=abs(dyqv(ij)) 460 538 ENDDO 539 461 540 462 541 c calcul des pentes aux points scalaires … … 475 554 IF (pole_nord) THEN 476 555 DO ij=1,iip1 477 dyq(ij,l)=qpns-q(ij+iip1,l )556 dyq(ij,l)=qpns-q(ij+iip1,l,iq) 478 557 ENDDO 479 558 … … 497 576 498 577 DO ij=1,iip1 499 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l )-qpsn578 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn 500 579 ENDDO 501 580 … … 633 712 DO ij=ijb,ije 634 713 IF(masse_adv_v(ij,l).gt.0) THEN 635 qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)* 636 , 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)) 714 qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)* 715 , 0.5*(1.-masse_adv_v(ij,l) 716 , /masse(ij+iip1,l,iq)) 637 717 ELSE 638 qbyv(ij,l)=q(ij,l )-dyq(ij,l)*639 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l ))718 qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)* 719 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) 640 720 ENDIF 641 721 qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l) … … 643 723 ENDDO 644 724 c$OMP END DO NOWAIT 725 726 ! CRisi: appel récursif de l'advection sur les fils. 727 ! Il faut faire ça avant d'avoir mis à jour q et masse 728 write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq) 729 730 ijb=ij_begin-2*iip1 731 ije=ij_end+2*iip1 732 if (pole_nord) ijb=ij_begin 733 if (pole_sud) ije=ij_end 734 735 if (nqfils(iq).gt.0) then 736 do ifils=1,nqdesc(iq) 737 iq2=iqfils(ifils,iq) 738 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 739 DO l=1,llm 740 DO ij=ijb,ije 741 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 742 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 743 enddo 744 enddo 745 c$OMP END DO NOWAIT 746 enddo !do ifils=1,nqdesc(iq) 747 748 do ifils=1,nqfils(iq) 749 iq2=iqfils(ifils,iq) 750 call vly_loc(Ratio,pente_max,masseq,qbyv,iq2) 751 enddo !do ifils=1,nqfils(iq) 752 endif !if (nqfils(iq).gt.0) then 753 ! end CRisi 645 754 646 755 ijb=ij_begin … … 652 761 DO l=1,llm 653 762 DO ij=ijb,ije 654 newmasse=masse(ij,l) 655 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l) 656 657 q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l)) 658 & /newmasse 659 masse(ij,l)=newmasse 660 ENDDO 763 newmasse=masse(ij,l,iq) 764 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l) 765 766 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l) 767 & -qbyv(ij-iip1,l))/newmasse 768 769 masse(ij,l,iq)=newmasse 770 771 ENDDO 772 773 661 774 c.-. ancienne version 662 775 c convpn=SSUM(iim,qbyv(1,l),1)/apoln … … 665 778 convpn=SSUM(iim,qbyv(1,l),1) 666 779 convmpn=ssum(iim,masse_adv_v(1,l),1) 667 massepn=ssum(iim,masse(1,l ),1)780 massepn=ssum(iim,masse(1,l,iq),1) 668 781 qpn=0. 669 782 do ij=1,iim 670 qpn=qpn+masse(ij,l )*q(ij,l)783 qpn=qpn+masse(ij,l,iq)*q(ij,l,iq) 671 784 enddo 672 785 qpn=(qpn+convpn)/(massepn+convmpn) 673 786 do ij=1,iip1 674 q(ij,l )=qpn787 q(ij,l,iq)=qpn 675 788 enddo 676 789 endif … … 683 796 convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) 684 797 convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1) 685 masseps=ssum(iim, masse(ip1jm+1,l ),1)798 masseps=ssum(iim, masse(ip1jm+1,l,iq),1) 686 799 qps=0. 687 800 do ij = ip1jm+1,ip1jmp1-1 688 qps=qps+masse(ij,l )*q(ij,l)801 qps=qps+masse(ij,l,iq)*q(ij,l,iq) 689 802 enddo 690 803 qps=(qps+convps)/(masseps+convmps) 691 804 do ij=ip1jm+1,ip1jmp1 692 q(ij,l )=qps805 q(ij,l,iq)=qps 693 806 enddo 694 807 endif … … 704 817 c DO ij = 1,iip1 705 818 c q(ij,l)=newq 706 c masse(ij,l )=newmasse*aire(ij)819 c masse(ij,l,iq)=newmasse*aire(ij) 707 820 c ENDDO 708 821 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) … … 714 827 c DO ij = ip1jm+1,ip1jmp1 715 828 c q(ij,l)=newq 716 c masse(ij,l )=newmasse*aire(ij)829 c masse(ij,l,iq)=newmasse*aire(ij) 717 830 c ENDDO 718 831 c._. fin nouvelle version 719 832 ENDDO 720 833 c$OMP END DO NOWAIT 834 835 ! retablir les fils en rapport de melange par rapport a l'air: 836 ijb=ij_begin 837 ije=ij_end 838 ! if (pole_nord) ijb=ij_begin 839 ! if (pole_sud) ije=ij_end 840 841 if (nqfils(iq).gt.0) then 842 do ifils=1,nqdesc(iq) 843 iq2=iqfils(ifils,iq) 844 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 845 DO l=1,llm 846 DO ij=ijb,ije 847 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 848 enddo 849 enddo 850 c$OMP END DO NOWAIT 851 enddo !do ifils=1,nqdesc(iq) 852 endif !if (nqfils(iq).gt.0) then 853 721 854 722 855 RETURN … … 725 858 726 859 727 SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x )860 SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x,iq) 728 861 c 729 862 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 739 872 USE parallel_lmdz 740 873 USE vlz_mod 874 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 741 875 IMPLICIT NONE 742 876 c … … 750 884 c Arguments: 751 885 c ---------- 752 REAL masse(ijb_u:ije_u,llm),pente_max 753 REAL q(ijb_u:ije_u,llm) 754 REAL w(ijb_u:ije_u,llm+1) 886 REAL masse(ijb_u:ije_u,llm,nqtot),pente_max 887 REAL q(ijb_u:ije_u,llm,nqtot) 888 REAL w(ijb_u:ije_u,llm+1,nqtot) 889 INTEGER iq 755 890 c 756 891 c Local … … 779 914 LOGICAL,SAVE :: first=.TRUE. 780 915 !$OMP THREADPRIVATE(first) 781 916 917 !REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi 918 ! Ces varibles doivent être déclarées en pointer et en save dans 919 ! vlz_loc si on veut qu'elles soient vues par tous les threads. 920 INTEGER ifils,iq2 ! CRisi 782 921 783 922 IF (first) THEN … … 787 926 c sens de W 788 927 928 write(*,*) 'vlsplt 926: entree dans vlz_loc, iq=',iq 789 929 #ifdef BIDON 790 930 IF(testcpu) THEN … … 799 939 DO l=2,llm 800 940 DO ij=ijb,ije 801 dzqw(ij,l)=q(ij,l-1 )-q(ij,l)941 dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq) 802 942 adzqw(ij,l)=abs(dzqw(ij,l)) 803 943 ENDDO … … 842 982 c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq 843 983 984 write(*,*) 'vlz 982,ijb,ije=',ijb,ije 844 985 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 845 986 DO l = 1,llm-1 846 987 do ij = ijb,ije 847 IF(w(ij,l+1).gt.0.) THEN 848 sigw=w(ij,l+1)/masse(ij,l+1) 849 wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1)) 988 IF(w(ij,l+1,iq).gt.0.) THEN 989 sigw=w(ij,l+1,iq)/masse(ij,l+1,iq) 990 wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l+1,iq) 991 : +0.5*(1.-sigw)*dzq(ij,l+1)) 850 992 ELSE 851 sigw=w(ij,l+1)/masse(ij,l) 852 wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l)) 993 sigw=w(ij,l+1,iq)/masse(ij,l,iq) 994 wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l,iq) 995 : -0.5*(1.+sigw)*dzq(ij,l)) 853 996 ENDIF 854 997 ENDDO 855 998 ENDDO 856 c$OMP END DO NOWAIT 999 c$OMP END DO NOWAIT 1000 write(*,*) 'vlz 1001' 857 1001 858 1002 c$OMP MASTER 859 1003 DO ij=ijb,ije 860 wq(ij,llm+1 )=0.861 wq(ij,1 )=0.1004 wq(ij,llm+1,iq)=0. 1005 wq(ij,1,iq)=0. 862 1006 ENDDO 863 1007 c$OMP END MASTER 864 1008 c$OMP BARRIER 865 1009 1010 ! CRisi: appel récursif de l'advection sur les fils. 1011 ! Il faut faire ça avant d'avoir mis à jour q et masse 1012 write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq) 1013 if (nqfils(iq).gt.0) then 1014 do ifils=1,nqdesc(iq) 1015 iq2=iqfils(ifils,iq) 1016 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1017 DO l=1,llm 1018 DO ij=ijb,ije 1019 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 1020 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 1021 wq(ij,l,iq2)=wq(ij,l,iq) 1022 enddo 1023 enddo 1024 c$OMP END DO NOWAIT 1025 enddo !do ifils=1,nqdesc(iq) 1026 c$OMP BARRIER 1027 1028 do ifils=1,nqfils(iq) 1029 iq2=iqfils(ifils,iq) 1030 call vlz_loc(Ratio,pente_max,masseq,wq,ijb_x,ije_x,iq2) 1031 enddo !do ifils=1,nqfils(iq) 1032 endif !if (nqfils(iq).gt.0) then 1033 ! end CRisi 1034 1035 ! CRisi: On rajoute ici une barrière car on veut être sur que tous les 1036 ! wq soient synchronisés 1037 write(*,*) 'vlz 1032' 1038 c$OMP BARRIER 866 1039 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 867 1040 DO l=1,llm 868 1041 DO ij=ijb,ije 869 newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l) 870 q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l)) 1042 newmasse=masse(ij,l,iq)+w(ij,l+1,iq)-w(ij,l,iq) 1043 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq) 1044 & +wq(ij,l+1,iq)-wq(ij,l,iq)) 871 1045 & /newmasse 872 masse(ij,l)=newmasse 873 ENDDO 874 ENDDO 875 c$OMP END DO NOWAIT 1046 masse(ij,l,iq)=newmasse 1047 ENDDO 1048 ENDDO 1049 c$OMP END DO NOWAIT 1050 1051 ! retablir les fils en rapport de melange par rapport a l'air: 1052 if (nqfils(iq).gt.0) then 1053 do ifils=1,nqdesc(iq) 1054 iq2=iqfils(ifils,iq) 1055 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1056 DO l=1,llm 1057 DO ij=ijb,ije 1058 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 1059 enddo 1060 enddo 1061 c$OMP END DO NOWAIT 1062 enddo !do ifils=1,nqdesc(iq) 1063 endif !if (nqfils(iq).gt.0) then 876 1064 877 1065 -
LMDZ5/trunk/libf/dyn3dmem/vlspltgen_loc.F
r1907 r2270 27 27 USE Write_Field_loc 28 28 USE VAMPIR 29 USE infotrac, ONLY : nqtot 29 ! CRisi: on rajoute variables utiles d'infotrac 30 USE infotrac, ONLY : nqtot,nqperes,nqdesc,nqfils,iqfils, 31 & ok_iso_verif 30 32 USE vlspltgen_mod 31 33 IMPLICIT NONE … … 64 66 REAL ptarg,pdelarg,foeew,zdelta 65 67 REAL tempe(ijb_u:ije_u) 66 INTEGER ijb,ije,iq 68 INTEGER ijb,ije,iq,iq2,ifils 67 69 LOGICAL, SAVE :: firstcall=.TRUE. 68 70 !$OMP THREADPRIVATE(firstcall) … … 150 152 ije=ij_end 151 153 154 DO iq=1,nqtot 152 155 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 153 156 DO l=1,llm 154 157 DO ij=ijb,ije 155 mw(ij,l )=w(ij,l) * zzw158 mw(ij,l,iq)=w(ij,l) * zzw 156 159 ENDDO 157 160 ENDDO 158 161 c$OMP END DO NOWAIT 159 162 ENDDO 163 164 DO iq=1,nqtot 160 165 c$OMP MASTER 161 166 DO ij=ijb,ije 162 mw(ij,llm+1 )=0.167 mw(ij,llm+1,iq)=0. 163 168 ENDDO 164 169 c$OMP END MASTER 170 ENDDO 165 171 166 172 c CALL SCOPY(ijp1llm,q,1,zq,1) … … 170 176 ije=ij_end 171 177 172 DO iq=1,nqtot 178 DO iq=1,nqtot 173 179 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 174 180 DO l=1,llm … … 179 185 ENDDO 180 186 181 #ifdef DEBUG_IO 187 #ifdef DEBUG_IO 182 188 CALL WriteField_u('mu',mu) 183 189 CALL WriteField_v('mv',mv) … … 186 192 #endif 187 193 194 ! verif temporaire 195 ijb=ij_begin-2*iip1 196 ije=ij_end+2*iip1 197 if (pole_nord) ijb=ij_begin 198 if (pole_sud) ije=ij_end 199 if (ok_iso_verif) then 200 call check_isotopes(zq,ijb,ije,'vlspltgen_loc 191') 201 endif !if (ok_iso_verif) then 202 188 203 c$OMP BARRIER 189 DO iq=1,nqtot 190 204 ! DO iq=1,nqtot 205 DO iq=1,nqperes ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air 206 write(*,*) 'vlspltgen 192: iq,iadv=',iq,iadv(iq) 207 #ifdef DEBUG_IO 208 CALL WriteField_u('zq',zq(:,:,iq)) 209 CALL WriteField_u('zm',zm(:,:,iq)) 210 #endif 211 if(iadv(iq) == 0) then 212 213 cycle 214 215 else if (iadv(iq)==10) then 216 217 #ifdef _ADV_HALO 218 ! CRisi: on ajoute les nombres de fils et tableaux des fils 219 ! On suppose qu'on ne peut advecter les fils que par le schéma 10. 220 call vlx_loc(zq,pente_max,zm,mu, 221 & ij_begin,ij_begin+2*iip1-1,iq) 222 call vlx_loc(zq,pente_max,zm,mu, 223 & ij_end-2*iip1+1,ij_end,iq) 224 #else 225 call vlx_loc(zq,pente_max,zm,mu, 226 & ij_begin,ij_end,iq) 227 #endif 228 229 c$OMP MASTER 230 call VTb(VTHallo) 231 c$OMP END MASTER 232 call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1) 233 call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1) 234 ! CRisi 235 do ifils=1,nqdesc(iq) 236 iq2=iqfils(ifils,iq) 237 call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1) 238 call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1) 239 enddo 240 241 c$OMP MASTER 242 call VTe(VTHallo) 243 c$OMP END MASTER 244 else if (iadv(iq)==14) then 245 246 #ifdef _ADV_HALO 247 call vlxqs_loc(zq,pente_max,zm,mu, 248 & qsat,ij_begin,ij_begin+2*iip1-1,iq) 249 call vlxqs_loc(zq,pente_max,zm,mu, 250 & qsat,ij_end-2*iip1+1,ij_end,iq) 251 #else 252 call vlxqs_loc(zq,pente_max,zm,mu, 253 & qsat,ij_begin,ij_end,iq) 254 #endif 255 256 c$OMP MASTER 257 call VTb(VTHallo) 258 c$OMP END MASTER 259 260 call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1) 261 call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1) 262 do ifils=1,nqdesc(iq) 263 iq2=iqfils(ifils,iq) 264 call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1) 265 call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1) 266 enddo 267 268 c$OMP MASTER 269 call VTe(VTHallo) 270 c$OMP END MASTER 271 else 272 273 stop 'vlspltgen_p : schema non parallelise' 274 275 endif 276 277 enddo !DO iq=1,nqperes 278 279 280 c$OMP BARRIER 281 c$OMP MASTER 282 call VTb(VTHallo) 283 c$OMP END MASTER 284 285 call SendRequest(MyRequest1) 286 287 c$OMP MASTER 288 call VTe(VTHallo) 289 c$OMP END MASTER 290 c$OMP BARRIER 291 292 ! verif temporaire 293 ijb=ij_begin-2*iip1 294 ije=ij_end+2*iip1 295 if (pole_nord) ijb=ij_begin 296 if (pole_sud) ije=ij_end 297 if (ok_iso_verif) then 298 call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 280') 299 endif !if (ok_iso_verif) then 300 301 do iq=1,nqperes 302 write(*,*) 'vlspltgen 279: iq=',iq 303 304 if(iadv(iq) == 0) then 305 306 cycle 307 308 else if (iadv(iq)==10) then 309 310 #ifdef _ADV_HALLO 311 call vlx_loc(zq,pente_max,zm,mu, 312 & ij_begin+2*iip1,ij_end-2*iip1,iq) 313 #endif 314 else if (iadv(iq)==14) then 315 #ifdef _ADV_HALLO 316 call vlxqs_loc(zq,pente_max,zm,mu, 317 & qsat,ij_begin+2*iip1,ij_end-2*iip1,iq) 318 #endif 319 else 320 321 stop 'vlspltgen_p : schema non parallelise' 322 323 endif 324 325 enddo 326 c$OMP BARRIER 327 c$OMP MASTER 328 call VTb(VTHallo) 329 c$OMP END MASTER 330 331 ! call WaitRecvRequest(MyRequest1) 332 ! call WaitSendRequest(MyRequest1) 333 c$OMP BARRIER 334 call WaitRequest(MyRequest1) 335 336 337 c$OMP MASTER 338 call VTe(VTHallo) 339 c$OMP END MASTER 340 c$OMP BARRIER 341 342 343 if (ok_iso_verif) then 344 call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 326') 345 endif !if (ok_iso_verif) then 346 if (ok_iso_verif) then 347 ijb=ij_begin-2*iip1 348 ije=ij_end+2*iip1 349 if (pole_nord) ijb=ij_begin 350 if (pole_sud) ije=ij_end 351 call check_isotopes(zq,ijb,ije,'vlspltgen_loc 336') 352 endif !if (ok_iso_verif) then 353 354 do iq=1,nqperes 355 write(*,*) 'vlspltgen 321: iq=',iq 191 356 #ifdef DEBUG_IO 192 357 CALL WriteField_u('zq',zq(:,:,iq)) 193 358 CALL WriteField_u('zm',zm(:,:,iq)) 194 359 #endif 360 195 361 if(iadv(iq) == 0) then 196 362 … … 198 364 199 365 else if (iadv(iq)==10) then 200 201 #ifdef _ADV_HALO 202 call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu, 203 & ij_begin,ij_begin+2*iip1-1) 204 call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu, 205 & ij_end-2*iip1+1,ij_end) 206 #else 207 call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu, 208 & ij_begin,ij_end) 209 #endif 210 211 c$OMP MASTER 212 call VTb(VTHallo) 213 c$OMP END MASTER 214 call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1) 215 call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1) 216 217 c$OMP MASTER 218 call VTe(VTHallo) 219 c$OMP END MASTER 366 367 call vly_loc(zq,pente_max,zm,mv,iq) 368 220 369 else if (iadv(iq)==14) then 221 222 #ifdef _ADV_HALO 223 call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu, 224 & qsat,ij_begin,ij_begin+2*iip1-1) 225 call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu, 226 & qsat,ij_end-2*iip1+1,ij_end) 227 #else 228 229 call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu, 230 & qsat,ij_begin,ij_end) 231 #endif 232 233 c$OMP MASTER 234 call VTb(VTHallo) 235 c$OMP END MASTER 236 237 call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1) 238 call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1) 239 240 c$OMP MASTER 241 call VTe(VTHallo) 242 c$OMP END MASTER 370 371 call vlyqs_loc(zq,pente_max,zm,mv, 372 & qsat,iq) 373 243 374 else 244 375 … … 246 377 247 378 endif 248 249 enddo 250 251 252 c$OMP BARRIER 253 c$OMP MASTER 254 call VTb(VTHallo) 255 c$OMP END MASTER 256 257 call SendRequest(MyRequest1) 258 259 c$OMP MASTER 260 call VTe(VTHallo) 261 c$OMP END MASTER 262 c$OMP BARRIER 263 do iq=1,nqtot 264 265 if(iadv(iq) == 0) then 266 267 cycle 268 269 else if (iadv(iq)==10) then 270 271 #ifdef _ADV_HALLO 272 call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu, 273 & ij_begin+2*iip1,ij_end-2*iip1) 274 #endif 275 else if (iadv(iq)==14) then 276 #ifdef _ADV_HALLO 277 call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu, 278 & qsat,ij_begin+2*iip1,ij_end-2*iip1) 279 #endif 280 else 281 282 stop 'vlspltgen_p : schema non parallelise' 283 284 endif 285 286 enddo 287 c$OMP BARRIER 288 c$OMP MASTER 289 call VTb(VTHallo) 290 c$OMP END MASTER 291 292 ! call WaitRecvRequest(MyRequest1) 293 ! call WaitSendRequest(MyRequest1) 294 c$OMP BARRIER 295 call WaitRequest(MyRequest1) 296 297 298 c$OMP MASTER 299 call VTe(VTHallo) 300 c$OMP END MASTER 301 c$OMP BARRIER 302 303 304 do iq=1,nqtot 379 380 enddo 381 382 if (ok_iso_verif) then 383 call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 357') 384 endif !if (ok_iso_verif) then 385 386 do iq=1,nqperes 387 write(*,*) 'vlspltgen 349: iq=',iq 305 388 #ifdef DEBUG_IO 306 389 CALL WriteField_u('zq',zq(:,:,iq)) 307 390 CALL WriteField_u('zm',zm(:,:,iq)) 308 391 #endif 392 if(iadv(iq) == 0) then 393 394 cycle 395 396 else if (iadv(iq)==10 .or. iadv(iq)==14 ) then 397 398 c$OMP BARRIER 399 #ifdef _ADV_HALLO 400 call vlz_loc(zq,pente_max,zm,mw, 401 & ij_begin,ij_begin+2*iip1-1,iq) 402 call vlz_loc(zq,pente_max,zm,mw, 403 & ij_end-2*iip1+1,ij_end,iq) 404 #else 405 call vlz_loc(zq,pente_max,zm,mw, 406 & ij_begin,ij_end,iq) 407 #endif 408 c$OMP BARRIER 409 410 c$OMP MASTER 411 call VTb(VTHallo) 412 c$OMP END MASTER 413 414 call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2) 415 call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2) 416 ! CRisi 417 do ifils=1,nqdesc(iq) 418 iq2=iqfils(ifils,iq) 419 call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest2) 420 call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest2) 421 enddo 422 c$OMP MASTER 423 call VTe(VTHallo) 424 c$OMP END MASTER 425 c$OMP BARRIER 426 else 427 428 stop 'vlspltgen_p : schema non parallelise' 429 430 endif 431 432 enddo 433 c$OMP BARRIER 434 435 c$OMP MASTER 436 call VTb(VTHallo) 437 c$OMP END MASTER 438 439 call SendRequest(MyRequest2) 440 441 c$OMP MASTER 442 call VTe(VTHallo) 443 c$OMP END MASTER 444 445 446 if (ok_iso_verif) then 447 call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 429') 448 endif !if (ok_iso_verif) then 449 450 c$OMP BARRIER 451 do iq=1,nqperes 452 write(*,*) 'vlspltgen 409: iq=',iq 453 309 454 if(iadv(iq) == 0) then 310 455 311 456 cycle 312 457 313 else if (iadv(iq)==10) then 314 315 call vly_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv) 316 317 else if (iadv(iq)==14) then 318 319 call vlyqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv, 320 & qsat) 321 458 else if (iadv(iq)==10 .or. iadv(iq)==14 ) then 459 c$OMP BARRIER 460 461 write(*,*) 'vlspltgen_loc 461' 462 #ifdef _ADV_HALLO 463 write(*,*) 'vlspltgen_loc 462' 464 call vlz_loc(zq,pente_max,zm,mw, 465 & ij_begin+2*iip1,ij_end-2*iip1,iq) 466 #endif 467 468 c$OMP BARRIER 322 469 else 323 470 … … 325 472 326 473 endif 327 328 enddo 329 330 331 do iq=1,nqtot 474 475 enddo 476 write(*,*) 'vlspltgen_loc 476' 477 478 c$OMP BARRIER 479 c$OMP MASTER 480 call VTb(VTHallo) 481 c$OMP END MASTER 482 483 ! call WaitRecvRequest(MyRequest2) 484 ! call WaitSendRequest(MyRequest2) 485 c$OMP BARRIER 486 CALL WaitRequest(MyRequest2) 487 488 c$OMP MASTER 489 call VTe(VTHallo) 490 c$OMP END MASTER 491 c$OMP BARRIER 492 493 494 write(*,*) 'vlspltgen_loc 494' 495 if (ok_iso_verif) then 496 call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 461') 497 endif !if (ok_iso_verif) then 498 499 do iq=1,nqperes 500 write(*,*) 'vlspltgen 449: iq=',iq 332 501 #ifdef DEBUG_IO 333 502 CALL WriteField_u('zq',zq(:,:,iq)) 334 503 CALL WriteField_u('zm',zm(:,:,iq)) 335 504 #endif 505 if(iadv(iq) == 0) then 506 507 cycle 508 509 else if (iadv(iq)==10) then 510 511 call vly_loc(zq,pente_max,zm,mv,iq) 512 513 else if (iadv(iq)==14) then 514 515 call vlyqs_loc(zq,pente_max,zm,mv, 516 & qsat,iq) 517 518 else 519 520 stop 'vlspltgen_p : schema non parallelise' 521 522 endif 523 524 enddo !do iq=1,nqperes 525 526 if (ok_iso_verif) then 527 call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 493') 528 endif !if (ok_iso_verif) then 529 530 do iq=1,nqperes 531 write(*,*) 'vlspltgen 477: iq=',iq 532 #ifdef DEBUG_IO 533 CALL WriteField_u('zq',zq(:,:,iq)) 534 CALL WriteField_u('zm',zm(:,:,iq)) 535 #endif 336 536 if(iadv(iq) == 0) then 337 537 338 538 cycle 339 539 340 else if (iadv(iq)==10 .or. iadv(iq)==14 ) then341 342 c$OMP BARRIER343 #ifdef _ADV_HALLO344 call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,345 & ij_begin,ij_begin+2*iip1-1)346 call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,347 & ij_end-2*iip1+1,ij_end)348 #else349 call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,350 & ij_begin,ij_end)351 #endif352 c$OMP BARRIER353 354 c$OMP MASTER355 call VTb(VTHallo)356 c$OMP END MASTER357 358 call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2)359 call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2)360 361 c$OMP MASTER362 call VTe(VTHallo)363 c$OMP END MASTER364 c$OMP BARRIER365 else366 367 stop 'vlspltgen_p : schema non parallelise'368 369 endif370 371 enddo372 c$OMP BARRIER373 374 c$OMP MASTER375 call VTb(VTHallo)376 c$OMP END MASTER377 378 call SendRequest(MyRequest2)379 380 c$OMP MASTER381 call VTe(VTHallo)382 c$OMP END MASTER383 384 c$OMP BARRIER385 do iq=1,nqtot386 387 if(iadv(iq) == 0) then388 389 cycle390 391 else if (iadv(iq)==10 .or. iadv(iq)==14 ) then392 c$OMP BARRIER393 394 #ifdef _ADV_HALLO395 call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,396 & ij_begin+2*iip1,ij_end-2*iip1)397 #endif398 399 c$OMP BARRIER400 else401 402 stop 'vlspltgen_p : schema non parallelise'403 404 endif405 406 enddo407 408 c$OMP BARRIER409 c$OMP MASTER410 call VTb(VTHallo)411 c$OMP END MASTER412 413 ! call WaitRecvRequest(MyRequest2)414 ! call WaitSendRequest(MyRequest2)415 c$OMP BARRIER416 CALL WaitRequest(MyRequest2)417 418 c$OMP MASTER419 call VTe(VTHallo)420 c$OMP END MASTER421 c$OMP BARRIER422 423 424 do iq=1,nqtot425 #ifdef DEBUG_IO426 CALL WriteField_u('zq',zq(:,:,iq))427 CALL WriteField_u('zm',zm(:,:,iq))428 #endif429 if(iadv(iq) == 0) then430 431 cycle432 433 540 else if (iadv(iq)==10) then 434 541 435 call vly_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv) 542 call vlx_loc(zq,pente_max,zm,mu, 543 & ij_begin,ij_end,iq) 436 544 437 545 else if (iadv(iq)==14) then 438 546 439 call vl yqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv,440 & qsat)547 call vlxqs_loc(zq,pente_max,zm,mu, 548 & qsat, ij_begin,ij_end,iq) 441 549 442 550 else 443 551 444 552 stop 'vlspltgen_p : schema non parallelise' 445 553 446 554 endif 447 555 448 enddo 449 450 451 do iq=1,nqtot 452 #ifdef DEBUG_IO 453 CALL WriteField_u('zq',zq(:,:,iq)) 454 CALL WriteField_u('zm',zm(:,:,iq)) 455 #endif 456 if(iadv(iq) == 0) then 457 458 cycle 459 460 else if (iadv(iq)==10) then 461 462 call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu, 463 & ij_begin,ij_end) 464 465 else if (iadv(iq)==14) then 466 467 call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu, 468 & qsat, ij_begin,ij_end) 469 470 else 471 472 stop 'vlspltgen_p : schema non parallelise' 473 474 endif 475 476 enddo 477 556 enddo !do iq=1,nqperes 557 558 write(*,*) 'vlspltgen 550: apres derniere serie de call vlx' 559 if (ok_iso_verif) then 560 call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 521') 561 endif !if (ok_iso_verif) then 478 562 479 563 ijb=ij_begin 480 564 ije=ij_end 565 write(*,*) 'vlspltgen_loc 557' 481 566 c$OMP BARRIER 482 567 483 568 write(*,*) 'vlspltgen_loc 559' 484 569 DO iq=1,nqtot 570 write(*,*) 'vlspltgen_loc 561, iq=',iq 485 571 #ifdef DEBUG_IO 486 572 CALL WriteField_u('zq',zq(:,:,iq)) … … 495 581 ENDDO 496 582 ENDDO 497 c$OMP END DO NOWAIT 583 c$OMP END DO NOWAIT 584 write(*,*) 'vlspltgen_loc 575' 498 585 499 586 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) … … 504 591 ENDDO 505 592 c$OMP END DO NOWAIT 506 507 ENDDO 593 write(*,*) 'vlspltgen_loc 583' 594 ENDDO !DO iq=1,nqtot 508 595 509 596 if (ok_iso_verif) then 597 call check_isotopes(q,ij_begin,ij_end,'vlspltgen_loc 557') 598 endif !if (ok_iso_verif) then 599 510 600 c$OMP BARRIER 511 601 … … 516 606 cc$OMP BARRIER 517 607 608 write(*,*) 'vlspltgen 597: sortie' 518 609 RETURN 519 610 END -
LMDZ5/trunk/libf/dyn3dmem/vlspltgen_mod.F90
r1907 r2270 2 2 3 3 REAL,POINTER,SAVE :: qsat(:,:) 4 REAL,POINTER,SAVE :: mu(:,:) 4 REAL,POINTER,SAVE :: mu(:,:) ! CRisi: on ajoute une dimension 5 5 REAL,POINTER,SAVE :: mv(:,:) 6 REAL,POINTER,SAVE :: mw(:,: )6 REAL,POINTER,SAVE :: mw(:,:,:) 7 7 REAL,POINTER,SAVE :: zm(:,:,:) 8 8 REAL,POINTER,SAVE :: zq(:,:,:) … … 25 25 CALL allocate_u(mu,llm,d) 26 26 CALL allocate_v(mv,llm,d) 27 CALL allocate_u(mw,llm+1, d)27 CALL allocate_u(mw,llm+1,nqtot,d) 28 28 CALL allocate_u(zm,llm,nqtot,d) 29 29 CALL allocate_u(zq,llm,nqtot,d) -
LMDZ5/trunk/libf/dyn3dmem/vlspltqs_loc.F
r1907 r2270 1 SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x )1 SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x,iq) 2 2 c 3 3 c Auteurs: P.Le Van, F.Hourdin, F.Forget 4 4 c 5 5 c ******************************************************************** 6 c Shema d' advection " pseudo amont " .6 c Shema d''advection " pseudo amont " . 7 7 c ******************************************************************** 8 8 c 9 9 c -------------------------------------------------------------------- 10 10 USE parallel_lmdz 11 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 11 12 IMPLICIT NONE 12 13 c … … 20 21 c Arguments: 21 22 c ---------- 22 REAL masse(ijb_u:ije_u,llm ),pente_max23 REAL masse(ijb_u:ije_u,llm,nqtot),pente_max 23 24 REAL u_m( ijb_u:ije_u,llm ) 24 REAL q(ijb_u:ije_u,llm )25 REAL q(ijb_u:ije_u,llm,nqtot) 25 26 REAL qsat(ijb_u:ije_u,llm) 27 INTEGER iq ! CRisi 26 28 c 27 29 c Local … … 36 38 REAL adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm) 37 39 REAL u_mq(ijb_u:ije_u,llm) 40 REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi 41 INTEGER ifils,iq2 ! CRisi 42 38 43 39 44 REAL SSUM … … 42 47 INTEGER ijb,ije,ijb_x,ije_x 43 48 49 write(*,*) 'vlspltqs 58: entree vlxqs_loc, iq,ijb_x=', 50 & iq,ijb_x 44 51 45 52 c calcul de la pente a droite et a gauche de la maille … … 65 72 DO l = 1, llm 66 73 DO ij=ijb,ije-1 67 dxqu(ij)=q(ij+1,l )-q(ij,l)74 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 68 75 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0' 69 c sigu(ij)=u_m(ij,l)/masse(ij,l )76 c sigu(ij)=u_m(ij,l)/masse(ij,l,iq) 70 77 ENDDO 71 78 DO ij=ijb+iip1-1,ije,iip1 … … 120 127 DO l = 1, llm 121 128 DO ij=ijb,ije-1 122 dxqu(ij)=q(ij+1,l )-q(ij,l)129 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 123 130 ENDDO 124 131 DO ij=ijb+iip1-1,ije,iip1 … … 179 186 DO l=1,llm 180 187 DO ij=ijb,ije-1 181 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l ),182 , 1.+u_m(ij,l)/masse(ij+1,l ),188 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq), 189 , 1.+u_m(ij,l)/masse(ij+1,l,iq), 183 190 , u_m(ij,l)) 184 191 zdum(ij,l)=0.5*zdum(ij,l) 185 192 u_mq(ij,l)=cvmgp( 186 , q(ij,l )+zdum(ij,l)*dxq(ij,l),187 , q(ij+1,l )-zdum(ij,l)*dxq(ij+1,l),193 , q(ij,l,iq)+zdum(ij,l)*dxq(ij,l), 194 , q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l), 188 195 , u_m(ij,l)) 189 196 u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l) … … 195 202 c on cumule le flux correspondant a toutes les mailles dont la masse 196 203 c au travers de la paroi pENDant le pas de temps. 197 c le rapport de melange de l' air advecte est min(q_vanleer, Qsat_downwind)204 c le rapport de melange de l''air advecte est min(q_vanleer, Qsat_downwind) 198 205 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 199 206 DO l=1,llm 200 207 DO ij=ijb,ije-1 201 208 IF (u_m(ij,l).gt.0.) THEN 202 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l )209 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq) 203 210 u_mq(ij,l)=u_m(ij,l)* 204 $ min(q(ij,l )+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))211 $ min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l)) 205 212 ELSE 206 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l )213 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq) 207 214 u_mq(ij,l)=u_m(ij,l)* 208 $ min(q(ij+1,l )-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))215 $ min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l)) 209 216 ENDIF 210 217 ENDDO … … 273 280 ENDDO 274 281 niju=iju 275 c PRINT*,'niju,nl',niju,nl(l)282 PRINT*,'vlxqs 280: niju,nl',niju,nl(l) 276 283 277 284 c traitement des mailles … … 285 292 i=ijq-(j-1)*iip1 286 293 c accumulation pour les mailles completements advectees 287 do while(zu_m.gt.masse(ijq,l)) 288 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l) 289 zu_m=zu_m-masse(ijq,l) 294 do while(zu_m.gt.masse(ijq,l,iq)) 295 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq) 296 & *masse(ijq,l,iq) 297 zu_m=zu_m-masse(ijq,l,iq) 290 298 i=mod(i-2+iim,iim)+1 291 299 ijq=(j-1)*iip1+i 292 300 ENDDO 293 301 c ajout de la maille non completement advectee 294 u_mq(ij,l)=u_mq(ij,l)+zu_m* 295 & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))302 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq) 303 & +0.5*(1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 296 304 ELSE 297 305 ijq=ij+1 298 306 i=ijq-(j-1)*iip1 299 307 c accumulation pour les mailles completements advectees 300 do while(-zu_m.gt.masse(ijq,l)) 301 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l) 302 zu_m=zu_m+masse(ijq,l) 308 do while(-zu_m.gt.masse(ijq,l,iq)) 309 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq) 310 & *masse(ijq,l,iq) 311 zu_m=zu_m+masse(ijq,l,iq) 303 312 i=mod(i,iim)+1 304 313 ijq=(j-1)*iip1+i 305 314 ENDDO 306 315 c ajout de la maille non completement advectee 307 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l )-308 & 0.5*(1.+zu_m/masse(ijq,l ))*dxq(ijq,l))316 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)- 317 & 0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 309 318 ENDIF 310 319 ENDDO … … 325 334 c$OMP END DO NOWAIT 326 335 336 ! CRisi: appel récursif de l'advection sur les fils. 337 ! Il faut faire ça avant d'avoir mis à jour q et masse 338 write(*,*) 'vlspltqs 336: iq,ijb_x,nqfils(iq)=', 339 & iq,ijb_x,nqfils(iq) 340 341 if (nqfils(iq).gt.0) then 342 do ifils=1,nqdesc(iq) 343 iq2=iqfils(ifils,iq) 344 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 345 DO l=1,llm 346 DO ij=ijb,ije 347 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 348 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 349 enddo 350 enddo 351 c$OMP END DO NOWAIT 352 enddo !do ifils=1,nqfils(iq) 353 do ifils=1,nqfils(iq) 354 iq2=iqfils(ifils,iq) 355 write(*,*) 'vlxqs 349: on appelle vlx pour iq2=',iq2 356 call vlx_loc(Ratio,pente_max,masseq,u_mq,ijb_x,ije_x,iq2) 357 enddo !do ifils=1,nqfils(iq) 358 endif !if (nqfils(iq).gt.0) then 359 ! end CRisi 360 361 write(*,*) 'vlspltqs 360: iq,ijb_x=',iq,ijb_x 362 327 363 c calcul des tendances 328 364 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 329 365 DO l=1,llm 330 366 DO ij=ijb+1,ije 331 new_m=masse(ij,l )+u_m(ij-1,l)-u_m(ij,l)332 q(ij,l )=(q(ij,l)*masse(ij,l)+367 new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l) 368 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+ 333 369 & u_mq(ij-1,l)-u_mq(ij,l)) 334 370 & /new_m 335 masse(ij,l )=new_m336 ENDDO 337 c Modif Fred 22 03 96 correction d' un bug (les scopy ci-dessous)371 masse(ij,l,iq)=new_m 372 ENDDO 373 c Modif Fred 22 03 96 correction d''un bug (les scopy ci-dessous) 338 374 DO ij=ijb+iip1-1,ije,iip1 339 q(ij-iim,l)=q(ij,l) 340 masse(ij-iim,l)=masse(ij,l) 341 ENDDO 342 ENDDO 343 c$OMP END DO NOWAIT 375 q(ij-iim,l,iq)=q(ij,l,iq) 376 masse(ij-iim,l,iq)=masse(ij,l,iq) 377 ENDDO 378 ENDDO 379 c$OMP END DO NOWAIT 380 381 write(*,*) 'vlspltqs 380: iq,ijb_x=',iq,ijb_x 382 383 ! retablir les fils en rapport de melange par rapport a l'air: 384 if (nqfils(iq).gt.0) then 385 do ifils=1,nqdesc(iq) 386 iq2=iqfils(ifils,iq) 387 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 388 DO l=1,llm 389 DO ij=ijb+1,ije 390 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 391 enddo 392 DO ij=ijb+iip1-1,ije,iip1 393 q(ij-iim,l,iq2)=q(ij,l,iq2) 394 enddo ! DO ij=ijb+iip1-1,ije,iip1 395 enddo 396 c$OMP END DO NOWAIT 397 enddo !do ifils=1,nqdesc(iq) 398 endif !if (nqfils(iq).gt.0) then 399 400 write(*,*) 'vlspltqs 399: iq,ijb_x=',iq,ijb_x 401 344 402 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) 345 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1 ),iip1,masse(iip2,1),iip1)403 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1,iq),iip1,masse(iip2,1,iq),iip1) 346 404 347 405 348 406 RETURN 349 407 END 350 SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat )408 SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat,iq) 351 409 c 352 410 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 361 419 c -------------------------------------------------------------------- 362 420 USE parallel_lmdz 421 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 363 422 IMPLICIT NONE 364 423 c … … 373 432 c Arguments: 374 433 c ---------- 375 REAL masse(ijb_u:ije_u,llm ),pente_max434 REAL masse(ijb_u:ije_u,llm,nqtot),pente_max 376 435 REAL masse_adv_v( ijb_v:ije_v,llm) 377 REAL q(ijb_u:ije_u,llm )436 REAL q(ijb_u:ije_u,llm,nqtot) 378 437 REAL qsat(ijb_u:ije_u,llm) 438 INTEGER iq ! CRisi 379 439 c 380 440 c Local … … 386 446 REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v) 387 447 REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u) 388 REAL qbyv(ijb_v:ije_v,llm )448 REAL qbyv(ijb_v:ije_v,llm,nqtot) 389 449 390 450 REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs … … 402 462 c 403 463 c 464 REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi 465 INTEGER ifils,iq2 ! CRisi 466 404 467 REAL SSUM 405 468 … … 407 470 INTEGER ijb,ije 408 471 472 ijb=ij_begin-2*iip1 473 ije=ij_end+2*iip1 474 if (pole_nord) ijb=ij_begin 475 if (pole_sud) ije=ij_end 476 ij=3525 477 l=3 478 if ((ij.ge.ijb).and.(ij.le.ije)) then 479 write(*,*) 'vlyqs 480: ij,l,iq,ijb,q(ij,l,:)=', 480 & ij,l,iq,ijb,q(ij,l,:) 481 endif 482 409 483 IF(first) THEN 410 484 PRINT*,'Shema Amont nouveau appele dans Vanleer ' 485 PRINT*,'vlyqs_loc, iq=',iq 411 486 first=.false. 412 487 do i=2,iip1 … … 439 514 if (pole_nord) then 440 515 DO i = 1, iim 441 airescb(i) = aire(i+ iip1) * q(i+ iip1,l )516 airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq) 442 517 ENDDO 443 518 qpns = SSUM( iim, airescb ,1 ) / airej2 … … 446 521 if (pole_sud) then 447 522 DO i = 1, iim 448 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l )523 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq) 449 524 ENDDO 450 525 qpsn = SSUM( iim, airesch ,1 ) / airejjm … … 460 535 461 536 DO ij=ijb,ije 462 dyqv(ij)=q(ij,l )-q(ij+iip1,l)537 dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq) 463 538 adyqv(ij)=abs(dyqv(ij)) 464 539 ENDDO … … 482 557 c calcul des pentes aux poles 483 558 DO ij=1,iip1 484 dyq(ij,l)=qpns-q(ij+iip1,l )559 dyq(ij,l)=qpns-q(ij+iip1,l,iq) 485 560 ENDDO 486 561 … … 513 588 514 589 DO ij=1,iip1 515 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l )-qpsn590 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn 516 591 ENDDO 517 592 … … 636 711 DO ij=ijb,ije 637 712 IF( masse_adv_v(ij,l).GT.0. ) THEN 638 qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l ) + 639 , dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))) 713 qbyv(ij,l,iq)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq ) + 714 , dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l) 715 , /masse(ij+iip1,l,iq))) 640 716 ELSE 641 qbyv(ij,l )= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *642 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l )) )717 qbyv(ij,l,iq)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) * 718 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) ) 643 719 ENDIF 644 qbyv(ij,l ) = masse_adv_v(ij,l)*qbyv(ij,l)720 qbyv(ij,l,iq) = masse_adv_v(ij,l)*qbyv(ij,l,iq) 645 721 ENDDO 646 722 ENDDO 647 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(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq) 728 729 ijb=ij_begin-2*iip1 730 ije=ij_end+2*iip1 731 if (pole_nord) ijb=ij_begin 732 if (pole_sud) ije=ij_end 733 734 if (nqfils(iq).gt.0) then 735 do ifils=1,nqdesc(iq) 736 iq2=iqfils(ifils,iq) 737 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 738 DO l=1,llm 739 DO ij=ijb,ije 740 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 741 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 742 enddo 743 enddo 744 c$OMP END DO NOWAIT 745 enddo !do ifils=1,nqdesc(iq) 746 do ifils=1,nqfils(iq) 747 iq2=iqfils(ifils,iq) 748 call vly_loc(Ratio,pente_max,masseq,qbyv,iq2) 749 enddo !do ifils=1,nqfils(iq) 750 endif !if (nqfils(iq).gt.0) then 751 752 753 ! end CRisi 648 754 649 755 ijb=ij_begin … … 655 761 DO l=1,llm 656 762 DO ij=ijb,ije 657 newmasse=masse(ij,l )763 newmasse=masse(ij,l,iq) 658 764 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l) 659 q(ij,l )=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))660 & /newmasse661 masse(ij,l )=newmasse765 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l,iq) 766 & -qbyv(ij-iip1,l,iq))/newmasse 767 masse(ij,l,iq)=newmasse 662 768 ENDDO 663 769 c.-. ancienne version … … 665 771 IF (pole_nord) THEN 666 772 667 convpn=SSUM(iim,qbyv(1,l ),1)/apoln773 convpn=SSUM(iim,qbyv(1,l,iq),1)/apoln 668 774 convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln 669 775 DO ij = 1,iip1 670 newmasse=masse(ij,l )+convmpn*aire(ij)671 q(ij,l )=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/776 newmasse=masse(ij,l,iq)+convmpn*aire(ij) 777 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/ 672 778 & newmasse 673 masse(ij,l )=newmasse779 masse(ij,l,iq)=newmasse 674 780 ENDDO 675 781 … … 678 784 IF (pole_sud) THEN 679 785 680 convps = -SSUM(iim,qbyv(ip1jm-iim,l ),1)/apols786 convps = -SSUM(iim,qbyv(ip1jm-iim,l,iq),iq,1)/apols 681 787 convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols 682 788 DO ij = ip1jm+1,ip1jmp1 683 newmasse=masse(ij,l )+convmps*aire(ij)684 q(ij,l )=(q(ij,l)*masse(ij,l)+convps*aire(ij))/789 newmasse=masse(ij,l,iq)+convmps*aire(ij) 790 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/ 685 791 & newmasse 686 masse(ij,l )=newmasse792 masse(ij,l,iq)=newmasse 687 793 ENDDO 688 794 … … 691 797 692 798 c._. nouvelle version 693 c convpn=SSUM(iim,qbyv(1,l ),1)799 c convpn=SSUM(iim,qbyv(1,l,iq),1) 694 800 c convmpn=ssum(iim,masse_adv_v(1,l),1) 695 c oldmasse=ssum(iim,masse(1,l ),1)801 c oldmasse=ssum(iim,masse(1,l,iq),1) 696 802 c newmasse=oldmasse+convmpn 697 c newq=(q(1,l )*oldmasse+convpn)/newmasse803 c newq=(q(1,l,iq)*oldmasse+convpn)/newmasse 698 804 c newmasse=newmasse/apoln 699 805 c DO ij = 1,iip1 700 c q(ij,l )=newq701 c masse(ij,l )=newmasse*aire(ij)806 c q(ij,l,iq)=newq 807 c masse(ij,l,iq)=newmasse*aire(ij) 702 808 c ENDDO 703 c convps=-SSUM(iim,qbyv(ip1jm-iim,l ),1)809 c convps=-SSUM(iim,qbyv(ip1jm-iim,l,iq),1) 704 810 c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1) 705 c oldmasse=ssum(iim,masse(ip1jm-iim,l ),1)811 c oldmasse=ssum(iim,masse(ip1jm-iim,l,iq),1) 706 812 c newmasse=oldmasse+convmps 707 c newq=(q(ip1jmp1,l )*oldmasse+convps)/newmasse813 c newq=(q(ip1jmp1,l,iq)*oldmasse+convps)/newmasse 708 814 c newmasse=newmasse/apols 709 815 c DO ij = ip1jm+1,ip1jmp1 710 c q(ij,l )=newq711 c masse(ij,l )=newmasse*aire(ij)816 c q(ij,l,iq)=newq 817 c masse(ij,l,iq)=newmasse*aire(ij) 712 818 c ENDDO 713 819 c._. fin nouvelle version 714 820 ENDDO 715 821 c$OMP END DO NOWAIT 822 823 ! retablir les fils en rapport de melange par rapport a l'air: 824 ijb=ij_begin 825 ije=ij_end 826 ! if (pole_nord) ijb=ij_begin+iip1 827 ! if (pole_sud) ije=ij_end-iip1 828 829 if (nqfils(iq).gt.0) then 830 do ifils=1,nqdesc(iq) 831 iq2=iqfils(ifils,iq) 832 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 833 DO l=1,llm 834 DO ij=ijb,ije 835 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 836 enddo 837 enddo 838 c$OMP END DO NOWAIT 839 enddo !do ifils=1,nqdesc(iq) 840 endif !if (nqfils(iq).gt.0) then 841 842 716 843 RETURN 717 844 END -
LMDZ5/trunk/libf/dyn3dmem/vlz_mod.F90
r1907 r2270 1 1 MODULE vlz_mod 2 2 3 REAL,POINTER,SAVE :: wq(:,: )3 REAL,POINTER,SAVE :: wq(:,:,:) 4 4 REAL,POINTER,SAVE :: dzq(:,:) 5 5 REAL,POINTER,SAVE :: dzqw(:,:) 6 6 REAL,POINTER,SAVE :: adzqw(:,:) 7 ! CRisi: pour les traceurs: 8 REAL,POINTER,SAVE :: masseq(:,:,:) 9 REAL,POINTER,SAVE :: Ratio(:,:,:) 7 10 8 11 CONTAINS … … 18 21 19 22 d=>distrib_vanleer 20 CALL allocate_u(wq,llm+1, d)23 CALL allocate_u(wq,llm+1,nqtot,d) 21 24 CALL allocate_u(dzq,llm,d) 22 25 CALL allocate_u(dzqw,llm,d) 23 26 CALL allocate_u(adzqw,llm,d) 27 if (nqdesc_tot.gt.0) then 28 CALL allocate_u(masseq,llm,nqtot,d) 29 CALL allocate_u(Ratio,llm,nqtot,d) 30 endif !if (nqdesc_tot.gt.0) then 24 31 25 32 END SUBROUTINE vlz_allocate … … 29 36 USE bands 30 37 USE parallel_lmdz 38 USE infotrac 31 39 IMPLICIT NONE 32 40 TYPE(distrib),INTENT(IN) :: dist … … 36 44 CALL switch_u(dzqw,distrib_vanleer,dist) 37 45 CALL switch_u(adzqw,distrib_vanleer,dist) 46 ! CRisi: 47 if (nqdesc_tot.gt.0) then 48 CALL switch_u(masseq,distrib_vanleer,dist) 49 CALL switch_u(Ratio,distrib_vanleer,dist) 50 endif !if (nqdesc_tot.gt.0) then 38 51 39 52 END SUBROUTINE vlz_switch_vanleer
Note: See TracChangeset
for help on using the changeset viewer.
