Changeset 2270
- Timestamp:
- May 7, 2015, 5:45:04 PM (9 years ago)
- Location:
- LMDZ5/trunk
- Files:
-
- 3 added
- 25 edited
-
arch/arch-X64_ADA.fcm (modified) (1 diff)
-
arch/arch-X64_CURIE.fcm (modified) (1 diff)
-
arch/arch-gfortran.fcm (modified) (1 diff)
-
arch/arch-gfortran_CICLAD.fcm (modified) (1 diff)
-
libf/dyn3d/advtrac.F90 (modified) (5 diffs)
-
libf/dyn3d/caladvtrac.F (modified) (4 diffs)
-
libf/dyn3d/check_isotopes.F (added)
-
libf/dyn3d/dynetat0.F (modified) (1 diff)
-
libf/dyn3d/iniacademic.F90 (modified) (2 diffs)
-
libf/dyn3d/leapfrog.F (modified) (12 diffs)
-
libf/dyn3d/qminimum.F (modified) (7 diffs)
-
libf/dyn3d/vlsplt.F (modified) (36 diffs)
-
libf/dyn3d/vlspltqs.F (modified) (26 diffs)
-
libf/dyn3d_common/infotrac.F90 (modified) (11 diffs)
-
libf/dyn3d_common/iso_verif_dyn.F (added)
-
libf/dyn3dmem/advtrac_loc.F (modified) (4 diffs)
-
libf/dyn3dmem/caladvtrac_loc.F (modified) (3 diffs)
-
libf/dyn3dmem/check_isotopes_loc.F (added)
-
libf/dyn3dmem/dynetat0_loc.F (modified) (2 diffs)
-
libf/dyn3dmem/iniacademic_loc.F90 (modified) (3 diffs)
-
libf/dyn3dmem/integrd_loc.F (modified) (8 diffs)
-
libf/dyn3dmem/leapfrog_loc.F (modified) (22 diffs)
-
libf/dyn3dmem/qminimum_loc.F (modified) (9 diffs)
-
libf/dyn3dmem/vlsplt_loc.F (modified) (43 diffs)
-
libf/dyn3dmem/vlspltgen_loc.F (modified) (12 diffs)
-
libf/dyn3dmem/vlspltgen_mod.F90 (modified) (2 diffs)
-
libf/dyn3dmem/vlspltqs_loc.F (modified) (26 diffs)
-
libf/dyn3dmem/vlz_mod.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/arch/arch-X64_ADA.fcm
r2115 r2270 5 5 %FPP_FLAGS -P -traditional -I/smplocal/pub/FFTW/3.3.3_dyn/include/ 6 6 %FPP_DEF NC_DOUBLE FFT_FFTW 7 %BASE_FFLAGS -auto - mcmodel=large -integer-size 32 -real-size 64 -align all7 %BASE_FFLAGS -auto -recursive -mcmodel=large -integer-size 32 -real-size 64 -align all 8 8 %PROD_FFLAGS -O2 -ip -fp-model strict -axAVX,SSE4.2 9 9 %DEV_FFLAGS -p -g -O1 -fpe0 -traceback -
LMDZ5/trunk/arch/arch-X64_CURIE.fcm
r2115 r2270 5 5 %FPP_FLAGS -P -traditional 6 6 %FPP_DEF NC_DOUBLE FFT_MKL 7 #%BASE_FFLAGS -xHost -i4 -r8 -auto -align all -I$(MKL_INC_DIR) -I$(MKL_INC_DIR)/intel64/lp647 #%BASE_FFLAGS -recursive -xHost -i4 -r8 -auto -align all -I$(MKL_INC_DIR) -I$(MKL_INC_DIR)/intel64/lp64 8 8 %BASE_FFLAGS -i4 -r8 -auto -align all -I$(MKL_INC_DIR) -I$(MKL_INC_DIR)/intel64/lp64 -fp-model strict 9 9 %PROD_FFLAGS -O2 -
LMDZ5/trunk/arch/arch-gfortran.fcm
r1907 r2270 5 5 %FPP_FLAGS -P -traditional 6 6 %FPP_DEF NC_DOUBLE 7 %BASE_FFLAGS -c -fdefault-real-8 7 %BASE_FFLAGS -c -fdefault-real-8 -frecursive 8 8 %PROD_FFLAGS -O3 9 9 %DEV_FFLAGS -O -
LMDZ5/trunk/arch/arch-gfortran_CICLAD.fcm
r2132 r2270 5 5 %FPP_FLAGS -P -traditional 6 6 %FPP_DEF NC_DOUBLE 7 %BASE_FFLAGS -c -fdefault-real-8 -fcray-pointer 7 %BASE_FFLAGS -c -fdefault-real-8 -fcray-pointer -frecursive 8 8 %PROD_FFLAGS -O3 9 9 %DEV_FFLAGS -O -Wall -fbounds-check -
LMDZ5/trunk/libf/dyn3d/advtrac.F90
r2239 r2270 9 9 ! M.A Filiberti (04/2002) 10 10 ! 11 USE infotrac, ONLY: nqtot, iadv 11 USE infotrac, ONLY: nqtot, iadv,nqperes,ok_iso_verif 12 12 USE control_mod, ONLY: iapp_tracvl, day_step 13 13 … … 223 223 ! Appel des sous programmes d'advection 224 224 !----------------------------------------------------------- 225 do iq=1,nqtot 225 226 if (ok_iso_verif) then 227 call check_isotopes_seq(q,1,ip1jmp1,'advtrac 162') 228 endif !if (ok_iso_verif) then 229 230 do iq=1,nqperes 226 231 ! call clock(t_initial) 227 232 if(iadv(iq) == 0) cycle … … 230 235 ! ---------------------------------------------------------------- 231 236 if(iadv(iq).eq.10) THEN 232 call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr) 233 237 ! CRisi: on fait passer tout q pour avoir acces aux fils 238 239 write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:) 240 call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq) 241 234 242 ! ---------------------------------------------------------------- 235 243 ! Schema "pseudo amont" + test sur humidite specifique … … 238 246 else if(iadv(iq).eq.14) then 239 247 ! 240 CALL vlspltqs( q(1,1,1), 2., massem, wg , & 241 pbarug,pbarvg,dtvr,p,pk,teta ) 248 write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:) 249 CALL vlspltqs( q, 2., massem, wg , & 250 pbarug,pbarvg,dtvr,p,pk,teta,iq) 251 242 252 ! ---------------------------------------------------------------- 243 253 ! Schema de Frederic Hourdin … … 388 398 end DO 389 399 400 if (ok_iso_verif) then 401 call check_isotopes_seq(q,ip1jmp1,'advtrac 397') 402 endif !if (ok_iso_verif) then 390 403 391 404 !------------------------------------------------------------------ -
LMDZ5/trunk/libf/dyn3d/caladvtrac.F
r1907 r2270 53 53 if (planet_type.eq."earth") then 54 54 C initialisation 55 dq(:,:,1:2)=q(:,:,1:2) 55 ! CRisi: il faut gérer tous les traceurs si on veut pouvoir faire des 56 ! isotopes 57 ! dq(:,:,1:2)=q(:,:,1:2) 58 dq(:,:,1:nqtot)=q(:,:,1:nqtot) 56 59 57 60 c test des valeurs minmax … … 81 84 ENDDO 82 85 ENDDO 83 84 CALL qminimum( q, 2, finmasse )86 87 CALL qminimum( q, nqtot, finmasse ) 85 88 86 89 CALL SCOPY ( ip1jmp1*llm, masse, 1, finmasse, 1 ) … … 92 95 dtvrtrac = iapp_tracvl * dtvr 93 96 c 94 DO iq = 1 , 297 DO iq = 1 , nqtot 95 98 DO l = 1 , llm 96 99 DO ij = 1,ip1jmp1 … … 105 108 if (planet_type.eq."earth") then 106 109 ! Earth-specific treatment for the first 2 tracers (water) 107 dq(:,:,1: 2)=0.110 dq(:,:,1:nqtot)=0. 108 111 endif ! of if (planet_type.eq."earth") 109 112 ENDIF ! of IF( iapptrac.EQ.iapp_tracvl ) -
LMDZ5/trunk/libf/dyn3d/dynetat0.F
r1930 r2270 297 297 write(lunout,*)" Il est donc initialise a zero" 298 298 q(:,:,:,iq)=0. 299 300 ! CRisi: pour les isotopes, on peut faire init théorique 301 ! distill de Rayleigh très simplifiée 302 if (ok_isotopes) then 303 if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.0)) then 304 q(:,:,:,iq)=q(:,:,:,iqpere(iq)) 305 & *tnat(iso_num(iq)) 306 & *(q(:,:,:,iqpere(iq))/30.e-3) 307 & **(alpha_ideal(iso_num(iq))-1) 308 endif 309 if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.1)) then 310 q(:,:,:,iq)=q(:,:,:,iqiso(iso_indnum(iq), 311 & phase_num(iq))) 312 endif 313 endif !if (ok_isotopes) then 299 314 ELSE 300 315 ierr = NF90_GET_VAR(nid, nvarid, q(:,:,:,iq)) -
LMDZ5/trunk/libf/dyn3d/iniacademic.F90
r2087 r2270 5 5 6 6 USE filtreg_mod, ONLY: inifilr 7 USE infotrac , ONLY : nqtot7 USE infotrac 8 8 USE control_mod, ONLY: day_step,planet_type 9 9 #ifdef CPP_IOIPSL … … 262 262 if (i == 2) q(:,:,i)=1.e-15 263 263 if (i.gt.2) q(:,:,i)=0. 264 265 ! CRisi: init des isotopes 266 ! distill de Rayleigh très simplifiée 267 if (ok_isotopes) then 268 if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then 269 q(:,:,i)=q(:,:,iqpere(i)) & 270 & *tnat(iso_num(i)) & 271 & *(q(:,:,iqpere(i))/30.e-3) & 272 & **(alpha_ideal(iso_num(i))-1) 273 endif 274 if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then 275 q(:,:,i)=q(:,:,iqiso(iso_indnum(i),phase_num(i))) 276 endif 277 endif !if (ok_isotopes) then 278 264 279 enddo 265 280 else 266 281 q(:,:,:)=0 267 282 endif ! of if (planet_type=="earth") 283 284 if (ok_iso_verif) then 285 call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc') 286 endif !if (ok_iso_verif) then 268 287 269 288 ! add random perturbation to temperature -
LMDZ5/trunk/libf/dyn3d/leapfrog.F
r2239 r2270 11 11 use IOIPSL 12 12 #endif 13 USE infotrac, ONLY: nqtot 13 USE infotrac, ONLY: nqtot,ok_iso_verif 14 14 USE guide_mod, ONLY : guide_main 15 15 USE write_field, ONLY: writefield … … 235 235 jH_cur = jH_cur - int(jH_cur) 236 236 237 if (ok_iso_verif) then 238 call check_isotopes_seq(q,ip1jmp1,'leapfrog 321') 239 endif !if (ok_iso_verif) then 237 240 238 241 #ifdef CPP_IOIPSL … … 265 268 ! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 ) 266 269 ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 ) 270 271 if (ok_iso_verif) then 272 call check_isotopes_seq(q,ip1jmp1,'leapfrog 400') 273 endif !if (ok_iso_verif) then 267 274 268 275 2 CONTINUE ! Matsuno backward or leapfrog step begins here … … 305 312 endif 306 313 314 315 if (ok_iso_verif) then 316 call check_isotopes_seq(q,ip1jmp1,'leapfrog 589') 317 endif !if (ok_iso_verif) then 318 307 319 c----------------------------------------------------------------------- 308 320 c calcul des tendances dynamiques: … … 321 333 c calcul des tendances advection des traceurs (dont l'humidite) 322 334 c ------------------------------------------------------------- 335 336 if (ok_iso_verif) then 337 call check_isotopes_seq(q,ip1jmp1, 338 & 'leapfrog 686: avant caladvtrac') 339 endif !if (ok_iso_verif) then 323 340 324 341 IF( forward. OR . leapf ) THEN … … 327 344 * p, masse, dq, teta, 328 345 . flxw, pk) 346 329 347 330 348 IF (offline) THEN … … 346 364 c ---------------------------------- 347 365 348 349 CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 , 366 if (ok_iso_verif) then 367 write(*,*) 'leapfrog 720' 368 call check_isotopes_seq(q,ip1jmp1,'leapfrog 756') 369 endif !if (ok_iso_verif) then 370 371 CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 , 350 372 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ) 351 373 ! $ finvmaold ) 352 374 375 if (ok_iso_verif) then 376 write(*,*) 'leapfrog 724' 377 call check_isotopes_seq(q,ip1jmp1,'leapfrog 762') 378 endif !if (ok_iso_verif) then 353 379 354 380 c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd) … … 516 542 CALL massdair(p,masse) 517 543 544 if (ok_iso_verif) then 545 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196') 546 endif !if (ok_iso_verif) then 518 547 519 548 c----------------------------------------------------------------------- … … 600 629 c preparation du pas d'integration suivant ...... 601 630 631 if (ok_iso_verif) then 632 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509') 633 endif !if (ok_iso_verif) then 634 602 635 IF ( .NOT.purmats ) THEN 603 636 c ........................................................ … … 657 690 ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin)) 658 691 692 if (ok_iso_verif) then 693 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584') 694 endif !if (ok_iso_verif) then 695 659 696 c----------------------------------------------------------------------- 660 697 c ecriture de la bande histoire: … … 735 772 ELSE ! of IF (.not.purmats) 736 773 774 if (ok_iso_verif) then 775 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664') 776 endif !if (ok_iso_verif) then 777 737 778 c ........................................................ 738 779 c .............. schema matsuno ............... … … 757 798 758 799 ELSE ! of IF(forward) i.e. backward step 800 801 if (ok_iso_verif) then 802 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698') 803 endif !if (ok_iso_verif) then 759 804 760 805 IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN -
LMDZ5/trunk/libf/dyn3d/qminimum.F
r1907 r2270 2 2 ! $Header$ 3 3 ! 4 SUBROUTINE qminimum( q,nq ,deltap )4 SUBROUTINE qminimum( q,nqtot,deltap ) 5 5 6 USE infotrac, ONLY: ok_isotopes,ntraciso,iqiso,ok_iso_verif 6 7 IMPLICIT none 7 8 c … … 13 14 #include "comvert.h" 14 15 c 15 INTEGER nq 16 REAL q(ip1jmp1,llm,nq ), deltap(ip1jmp1,llm)16 INTEGER nqtot 17 REAL q(ip1jmp1,llm,nqtot), deltap(ip1jmp1,llm) 17 18 c 18 19 INTEGER iq_vap, iq_liq … … 30 31 INTEGER i, k, iq 31 32 REAL zx_defau, zx_abc, zx_pump(ip1jmp1), pompe 33 34 real zx_defau_diag(ip1jmp1,llm,2) 35 real q_follow(ip1jmp1,llm,2) 32 36 c 33 37 REAL SSUM … … 36 40 SAVE imprim 37 41 DATA imprim /0/ 42 !INTEGER ijb,ije 43 !INTEGER Index_pump(ij_end-ij_begin+1) 44 !INTEGER nb_pump 45 INTEGER ixt 38 46 c 39 47 c Quand l'eau liquide est trop petite (ou negative), on prend … … 41 49 c (sans changer la temperature !) 42 50 c 51 52 if (ok_iso_verif) then 53 call check_isotopes_seq(q,ip1jmp1,'qminimum 52') 54 endif !if (ok_iso_verif) then 55 56 zx_defau_diag(:,:,:)=0.0 57 q_follow(:,:,:)=q(:,:,:) 43 58 DO 1000 k = 1, llm 44 59 DO 1040 i = 1, ip1jmp1 45 60 if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then 61 62 if (ok_isotopes) then 63 zx_defau_diag(i,k,iq_liq)=AMAX1 64 : ( seuil_liq - q(i,k,iq_liq), 0.0 ) 65 endif !if (ok_isotopes) then 66 46 67 q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq 47 68 q(i,k,iq_liq) = seuil_liq … … 59 80 DO i = 1, ip1jmp1 60 81 if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then 82 83 if (ok_isotopes) then 84 zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 0.0 ) 85 endif !if (ok_isotopes) then 86 61 87 q(i,k-1,iq) = q(i,k-1,iq) - ( seuil_vap - q(i,k,iq) ) * 62 88 & deltap(i,k) / deltap(i,k-1) … … 83 109 ENDDO 84 110 ENDIF 111 112 write(*,*) 'qminimum 128' 113 if (ok_isotopes) then 114 ! CRisi: traiter de même les traceurs d'eau 115 ! Mais il faut les prendre à l'envers pour essayer de conserver la 116 ! masse. 117 ! 1) pompage dans le sol 118 ! On suppose que ce pompage se fait sans isotopes -> on ne modifie 119 ! rien ici et on croise les doigts pour que ça ne soit pas trop 120 ! génant 121 DO i = 1,ip1jmp1 122 if (zx_pump(i).gt.0.0) then 123 q_follow(i,1,iq_vap)=q_follow(i,1,iq_vap)+zx_pump(i) 124 endif !if (zx_pump(i).gt.0.0) then 125 enddo !DO i = 1,ip1jmp1 126 127 ! 2) transfert de vap vers les couches plus hautes 128 !write(*,*) 'qminimum 139' 129 do k=2,llm 130 DO i = 1,ip1jmp1 131 if (zx_defau_diag(i,k,iq_vap).gt.0.0) then 132 ! on ajoute la vapeur en k 133 do ixt=1,ntraciso 134 q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap)) 135 : +zx_defau_diag(i,k,iq_vap) 136 : *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap) 137 138 ! et on la retranche en k-1 139 q(i,k-1,iqiso(ixt,iq_vap))=q(i,k-1,iqiso(ixt,iq_vap)) 140 : -zx_defau_diag(i,k,iq_vap) 141 : *deltap(i,k)/deltap(i,k-1) 142 : *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap) 143 144 enddo !do ixt=1,niso 145 q_follow(i,k,iq_vap)= q_follow(i,k,iq_vap) 146 : +zx_defau_diag(i,k,iq_vap) 147 q_follow(i,k-1,iq_vap)= q_follow(i,k-1,iq_vap) 148 : -zx_defau_diag(i,k,iq_vap) 149 : *deltap(i,k)/deltap(i,k-1) 150 endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then 151 enddo !DO i = 1, ip1jmp1 152 enddo !do k=2,llm 153 154 if (ok_iso_verif) then 155 call check_isotopes_seq(q,ip1jmp1,'qminimum 168') 156 endif !if (ok_iso_verif) then 157 158 159 ! 3) transfert d'eau de la vapeur au liquide 160 !write(*,*) 'qminimum 164' 161 do k=1,llm 162 DO i = 1,ip1jmp1 163 if (zx_defau_diag(i,k,iq_liq).gt.0.0) then 164 165 ! on ajoute eau liquide en k en k 166 do ixt=1,ntraciso 167 q(i,k,iqiso(ixt,iq_liq))=q(i,k,iqiso(ixt,iq_liq)) 168 : +zx_defau_diag(i,k,iq_liq) 169 : *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap) 170 ! et on la retranche à la vapeur en k 171 q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap)) 172 : -zx_defau_diag(i,k,iq_liq) 173 : *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap) 174 enddo !do ixt=1,niso 175 q_follow(i,k,iq_liq)= q_follow(i,k,iq_liq) 176 : +zx_defau_diag(i,k,iq_liq) 177 q_follow(i,k,iq_vap)= q_follow(i,k,iq_vap) 178 : -zx_defau_diag(i,k,iq_liq) 179 endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then 180 enddo !DO i = 1, ip1jmp1 181 enddo !do k=2,llm 182 183 if (ok_iso_verif) then 184 call check_isotopes_seq(q,ip1jmp1,'qminimum 197') 185 endif !if (ok_iso_verif) then 186 187 endif !if (ok_isotopes) then 188 !write(*,*) 'qminimum 188' 189 85 190 c 86 191 RETURN -
LMDZ5/trunk/libf/dyn3d/vlsplt.F
r1907 r2270 3 3 c 4 4 5 SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt) 5 SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq) 6 USE infotrac, ONLY: nqtot,nqdesc,iqfils 6 7 c 7 8 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 32 33 c REAL masse(iip1,jjp1,llm),pente_max 33 34 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm) 34 REAL q(ip1jmp1,llm )35 REAL q(ip1jmp1,llm,nqtot) 35 36 c REAL q(iip1,jjp1,llm) 36 37 REAL w(ip1jmp1,llm),pdt 38 INTEGER iq ! CRisi 37 39 c 38 40 c Local … … 42 44 INTEGER ijlqmin,iqmin,jqmin,lqmin 43 45 c 44 REAL zm(ip1jmp1,llm ),newmasse46 REAL zm(ip1jmp1,llm,nqtot),newmasse 45 47 REAL mu(ip1jmp1,llm) 46 48 REAL mv(ip1jm,llm) 47 49 REAL mw(ip1jmp1,llm+1) 48 REAL zq(ip1jmp1,llm ),zz50 REAL zq(ip1jmp1,llm,nqtot),zz 49 51 REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm) 50 52 REAL second,temps0,temps1,temps2,temps3 … … 55 57 SAVE temps1,temps2,temps3 56 58 INTEGER iminn,imaxx 59 INTEGER ifils,iq2 ! CRisi 57 60 58 61 REAL qmin,qmax … … 79 82 mw(ij,llm+1)=0. 80 83 ENDDO 81 82 CALL SCOPY(ijp1llm,q,1,zq,1) 83 CALL SCOPY(ijp1llm,masse,1,zm,1) 84 85 CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1) 86 CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1) 87 88 if (nqdesc(iq).gt.0) then 89 do ifils=1,nqdesc(iq) 90 iq2=iqfils(ifils,iq) 91 CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1) 92 enddo 93 endif !if (nqfils(iq).gt.0) then 84 94 85 95 cprint*,'Entree vlx1' 86 96 c call minmaxq(zq,qmin,qmax,'avant vlx ') 87 call vlx(zq,pente_max,zm,mu )97 call vlx(zq,pente_max,zm,mu,iq) 88 98 cprint*,'Sortie vlx1' 89 99 c call minmaxq(zq,qmin,qmax,'apres vlx1 ') 90 100 91 101 c print*,'Entree vly1' 92 call vly(zq,pente_max,zm,mv) 102 103 call vly(zq,pente_max,zm,mv,iq) 93 104 c call minmaxq(zq,qmin,qmax,'apres vly1 ') 94 105 cprint*,'Sortie vly1' 95 call vlz(zq,pente_max,zm,mw )106 call vlz(zq,pente_max,zm,mw,iq) 96 107 c call minmaxq(zq,qmin,qmax,'apres vlz ') 97 108 98 109 99 call vly(zq,pente_max,zm,mv )110 call vly(zq,pente_max,zm,mv,iq) 100 111 c call minmaxq(zq,qmin,qmax,'apres vly ') 101 112 102 113 103 call vlx(zq,pente_max,zm,mu )114 call vlx(zq,pente_max,zm,mu,iq) 104 115 c call minmaxq(zq,qmin,qmax,'apres vlx2 ') 105 116 … … 107 118 DO l=1,llm 108 119 DO ij=1,ip1jmp1 109 q(ij,l )=zq(ij,l)120 q(ij,l,iq)=zq(ij,l,iq) 110 121 ENDDO 111 122 DO ij=1,ip1jm+1,iip1 112 q(ij+iim,l)=q(ij,l) 113 ENDDO 114 ENDDO 123 q(ij+iim,l,iq)=q(ij,l,iq) 124 ENDDO 125 ENDDO 126 ! CRisi: aussi pour les fils 127 if (nqdesc(iq).gt.0) then 128 do ifils=1,nqdesc(iq) 129 iq2=iqfils(ifils,iq) 130 DO l=1,llm 131 DO ij=1,ip1jmp1 132 q(ij,l,iq2)=zq(ij,l,iq2) 133 ENDDO 134 DO ij=1,ip1jm+1,iip1 135 q(ij+iim,l,iq2)=q(ij,l,iq2) 136 ENDDO 137 ENDDO 138 enddo !do ifils=1,nqdesc(iq) 139 endif ! if (nqdesc(iq).gt.0) then 115 140 116 141 RETURN 117 142 END 118 SUBROUTINE vlx(q,pente_max,masse,u_m) 143 SUBROUTINE vlx(q,pente_max,masse,u_m,iq) 144 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 119 145 120 146 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 139 165 c Arguments: 140 166 c ---------- 141 REAL masse(ip1jmp1,llm ),pente_max167 REAL masse(ip1jmp1,llm,nqtot),pente_max 142 168 REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm) 143 REAL q(ip1jmp1,llm )169 REAL q(ip1jmp1,llm,nqtot) 144 170 REAL w(ip1jmp1,llm) 171 INTEGER iq ! CRisi 145 172 c 146 173 c Local … … 155 182 REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm) 156 183 REAL u_mq(ip1jmp1,llm) 184 185 ! CRisi 186 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) 187 INTEGER ifils,iq2 ! CRisi 157 188 158 189 Logical extremum,first,testcpu … … 188 219 DO l = 1, llm 189 220 DO ij=iip2,ip1jm-1 190 dxqu(ij)=q(ij+1,l )-q(ij,l)221 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 191 222 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0' 192 c sigu(ij)=u_m(ij,l)/masse(ij,l )223 c sigu(ij)=u_m(ij,l)/masse(ij,l,iq) 193 224 ENDDO 194 225 DO ij=iip1+iip1,ip1jm,iip1 … … 243 274 DO l = 1, llm 244 275 DO ij=iip2,ip1jm-1 245 dxqu(ij)=q(ij+1,l )-q(ij,l)276 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 246 277 ENDDO 247 278 DO ij=iip1+iip1,ip1jm,iip1 … … 285 316 DO l=1,llm 286 317 DO ij=iip2,ip1jm-1 287 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l ),288 , 1.+u_m(ij,l)/masse(ij+1,l ),318 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq), 319 , 1.+u_m(ij,l)/masse(ij+1,l,iq), 289 320 , u_m(ij,l)) 290 321 zdum(ij,l)=0.5*zdum(ij,l) 291 322 u_mq(ij,l)=cvmgp( 292 , q(ij,l )+zdum(ij,l)*dxq(ij,l),293 , q(ij+1,l )-zdum(ij,l)*dxq(ij+1,l),323 , q(ij,l,iq)+zdum(ij,l)*dxq(ij,l), 324 , q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l), 294 325 , u_m(ij,l)) 295 326 u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l) … … 303 334 DO l=1,llm 304 335 DO ij=iip2,ip1jm-1 305 c print*,'masse(',ij,')=',masse(ij,l )336 c print*,'masse(',ij,')=',masse(ij,l,iq) 306 337 IF (u_m(ij,l).gt.0.) THEN 307 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l )308 u_mq(ij,l)=u_m(ij,l)*(q(ij,l )+0.5*zdum(ij,l)*dxq(ij,l))338 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq) 339 u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l)) 309 340 ELSE 310 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l) 311 u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l)) 341 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq) 342 u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq) 343 & -0.5*zdum(ij,l)*dxq(ij+1,l)) 312 344 ENDIF 313 345 ENDDO … … 379 411 i=ijq-(j-1)*iip1 380 412 c accumulation pour les mailles completements advectees 381 do while(zu_m.gt.masse(ijq,l)) 382 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l) 383 zu_m=zu_m-masse(ijq,l) 413 do while(zu_m.gt.masse(ijq,l,iq)) 414 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq) 415 & *masse(ijq,l,iq) 416 zu_m=zu_m-masse(ijq,l,iq) 384 417 i=mod(i-2+iim,iim)+1 385 418 ijq=(j-1)*iip1+i … … 387 420 c ajout de la maille non completement advectee 388 421 u_mq(ij,l)=u_mq(ij,l)+zu_m* 389 & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l)) 422 & (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq)) 423 & *dxq(ijq,l)) 390 424 ELSE 391 425 ijq=ij+1 392 426 i=ijq-(j-1)*iip1 393 427 c accumulation pour les mailles completements advectees 394 do while(-zu_m.gt.masse(ijq,l)) 395 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l) 396 zu_m=zu_m+masse(ijq,l) 428 do while(-zu_m.gt.masse(ijq,l,iq)) 429 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq) 430 & *masse(ijq,l,iq) 431 zu_m=zu_m+masse(ijq,l,iq) 397 432 i=mod(i,iim)+1 398 433 ijq=(j-1)*iip1+i 399 434 ENDDO 400 435 c ajout de la maille non completement advectee 401 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l )-402 & 0.5*(1.+zu_m/masse(ijq,l ))*dxq(ijq,l))436 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)- 437 & 0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 403 438 ENDIF 404 439 ENDDO … … 417 452 ENDDO 418 453 454 ! CRisi: appel récursif de l'advection sur les fils. 455 ! Il faut faire ça avant d'avoir mis à jour q et masse 456 write(*,*) 'vlsplt 326: iq,nqfils(iq)=',iq,nqfils(iq) 457 458 if (nqdesc(iq).gt.0) then 459 do ifils=1,nqdesc(iq) 460 iq2=iqfils(ifils,iq) 461 DO l=1,llm 462 DO ij=iip2,ip1jm 463 ! On a besoin de q et masse seulement entre iip2 et ip1jm 464 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 465 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 466 enddo 467 enddo 468 enddo !do ifils=1,nqdesc(iq) 469 do ifils=1,nqfils(iq) 470 iq2=iqfils(ifils,iq) 471 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 472 enddo !do ifils=1,nqfils(iq) 473 endif !if (nqfils(iq).gt.0) then 474 ! end CRisi 475 419 476 420 477 c calcul des tENDances … … 422 479 DO l=1,llm 423 480 DO ij=iip2+1,ip1jm 424 new_m=masse(ij,l )+u_m(ij-1,l)-u_m(ij,l)425 q(ij,l )=(q(ij,l)*masse(ij,l)+481 new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l) 482 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+ 426 483 & u_mq(ij-1,l)-u_mq(ij,l)) 427 484 & /new_m 428 masse(ij,l )=new_m485 masse(ij,l,iq)=new_m 429 486 ENDDO 430 487 c ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous) 431 488 DO ij=iip1+iip1,ip1jm,iip1 432 q(ij-iim,l)=q(ij,l) 433 masse(ij-iim,l)=masse(ij,l) 434 ENDDO 435 ENDDO 489 q(ij-iim,l,iq)=q(ij,l,iq) 490 masse(ij-iim,l,iq)=masse(ij,l,iq) 491 ENDDO 492 ENDDO 493 494 ! retablir les fils en rapport de melange par rapport a l'air: 495 ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio 496 ! puis on boucle en longitude 497 if (nqdesc(iq).gt.0) then 498 do ifils=1,nqdesc(iq) 499 iq2=iqfils(ifils,iq) 500 DO l=1,llm 501 DO ij=iip2+1,ip1jm 502 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 503 enddo 504 DO ij=iip1+iip1,ip1jm,iip1 505 q(ij-iim,l,iq2)=q(ij,l,iq2) 506 enddo ! DO ij=ijb+iip1-1,ije,iip1 507 enddo !DO l=1,llm 508 enddo !do ifils=1,nqdesc(iq) 509 endif !if (nqfils(iq).gt.0) then 510 436 511 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) 437 512 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1) … … 440 515 RETURN 441 516 END 442 SUBROUTINE vly(q,pente_max,masse,masse_adv_v) 517 SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq) 518 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 443 519 c 444 520 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 464 540 c Arguments: 465 541 c ---------- 466 REAL masse(ip1jmp1,llm ),pente_max542 REAL masse(ip1jmp1,llm,nqtot),pente_max 467 543 REAL masse_adv_v( ip1jm,llm) 468 REAL q(ip1jmp1,llm), dq( ip1jmp1,llm) 544 REAL q(ip1jmp1,llm,nqtot), dq( ip1jmp1,llm) 545 INTEGER iq ! CRisi 469 546 c 470 547 c Local … … 491 568 SAVE sinlon,coslon,sinlondlon,coslondlon 492 569 SAVE airej2,airejjm 570 571 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi 572 INTEGER ifils,iq2 ! CRisi 573 493 574 c 494 575 c … … 497 578 DATA first,testcpu/.true.,.false./ 498 579 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./ 580 581 write(*,*) 'vly 578: entree, iq=',iq 499 582 500 583 IF(first) THEN … … 529 612 530 613 DO i = 1, iim 531 airescb(i) = aire(i+ iip1) * q(i+ iip1,l )532 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l )614 airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq) 615 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq) 533 616 ENDDO 534 617 qpns = SSUM( iim, airescb ,1 ) / airej2 … … 538 621 539 622 DO ij=1,ip1jm 540 dyqv(ij)=q(ij,l )-q(ij+iip1,l)623 dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq) 541 624 adyqv(ij)=abs(dyqv(ij)) 542 625 ENDDO … … 553 636 554 637 DO ij=1,iip1 555 dyq(ij,l)=qpns-q(ij+iip1,l )556 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l )-qpsn638 dyq(ij,l)=qpns-q(ij+iip1,l,iq) 639 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn 557 640 ENDDO 558 641 … … 675 758 ENDDO 676 759 760 write(*,*) 'vly 756' 677 761 DO l=1,llm 678 762 DO ij=1,ip1jm 679 763 IF(masse_adv_v(ij,l).gt.0) THEN 680 qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)* 681 , 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)) 764 qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)* 765 , 0.5*(1.-masse_adv_v(ij,l) 766 , /masse(ij+iip1,l,iq)) 682 767 ELSE 683 qbyv(ij,l)=q(ij,l)-dyq(ij,l)* 684 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) 768 qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)* 769 , 0.5*(1.+masse_adv_v(ij,l) 770 , /masse(ij,l,iq)) 685 771 ENDIF 686 772 qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l) … … 688 774 ENDDO 689 775 776 ! CRisi: appel récursif de l'advection sur les fils. 777 ! Il faut faire ça avant d'avoir mis à jour q et masse 778 write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq) 779 780 if (nqfils(iq).gt.0) then 781 do ifils=1,nqdesc(iq) 782 iq2=iqfils(ifils,iq) 783 DO l=1,llm 784 DO ij=1,ip1jmp1 785 ! attention, chaque fils doit avoir son masseq, sinon, le 1er 786 ! fils ecrase le masseq de ses freres. 787 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 788 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 789 enddo 790 enddo 791 enddo !do ifils=1,nqdesc(iq) 792 793 do ifils=1,nqfils(iq) 794 iq2=iqfils(ifils,iq) 795 call vly(Ratio,pente_max,masseq,qbyv,iq2) 796 enddo !do ifils=1,nqfils(iq) 797 endif !if (nqfils(iq).gt.0) then 690 798 691 799 DO l=1,llm 692 800 DO ij=iip2,ip1jm 693 newmasse=masse(ij,l )801 newmasse=masse(ij,l,iq) 694 802 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l) 695 q(ij,l )=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))696 & /newmasse697 masse(ij,l )=newmasse803 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l) 804 & -qbyv(ij-iip1,l))/newmasse 805 masse(ij,l,iq)=newmasse 698 806 ENDDO 699 807 c.-. ancienne version … … 703 811 convpn=SSUM(iim,qbyv(1,l),1) 704 812 convmpn=ssum(iim,masse_adv_v(1,l),1) 705 massepn=ssum(iim,masse(1,l ),1)813 massepn=ssum(iim,masse(1,l,iq),1) 706 814 qpn=0. 707 815 do ij=1,iim 708 qpn=qpn+masse(ij,l )*q(ij,l)816 qpn=qpn+masse(ij,l,iq)*q(ij,l,iq) 709 817 enddo 710 818 qpn=(qpn+convpn)/(massepn+convmpn) 711 819 do ij=1,iip1 712 q(ij,l )=qpn820 q(ij,l,iq)=qpn 713 821 enddo 714 822 … … 718 826 convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) 719 827 convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1) 720 masseps=ssum(iim, masse(ip1jm+1,l ),1)828 masseps=ssum(iim, masse(ip1jm+1,l,iq),1) 721 829 qps=0. 722 830 do ij = ip1jm+1,ip1jmp1-1 723 qps=qps+masse(ij,l )*q(ij,l)831 qps=qps+masse(ij,l,iq)*q(ij,l,iq) 724 832 enddo 725 833 qps=(qps+convps)/(masseps+convmps) 726 834 do ij=ip1jm+1,ip1jmp1 727 q(ij,l )=qps835 q(ij,l,iq)=qps 728 836 enddo 729 837 … … 739 847 c DO ij = 1,iip1 740 848 c q(ij,l)=newq 741 c masse(ij,l )=newmasse*aire(ij)849 c masse(ij,l,iq)=newmasse*aire(ij) 742 850 c ENDDO 743 851 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) … … 749 857 c DO ij = ip1jm+1,ip1jmp1 750 858 c q(ij,l)=newq 751 c masse(ij,l )=newmasse*aire(ij)859 c masse(ij,l,iq)=newmasse*aire(ij) 752 860 c ENDDO 753 861 c._. fin nouvelle version 754 862 ENDDO 863 864 ! retablir les fils en rapport de melange par rapport a l'air: 865 if (nqfils(iq).gt.0) then 866 do ifils=1,nqdesc(iq) 867 iq2=iqfils(ifils,iq) 868 DO l=1,llm 869 DO ij=1,ip1jmp1 870 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 871 enddo 872 enddo 873 enddo !do ifils=1,nqdesc(iq) 874 endif !if (nqfils(iq).gt.0) then 875 876 write(*,*) 'vly 853: sortie' 755 877 756 878 RETURN 757 879 END 758 SUBROUTINE vlz(q,pente_max,masse,w) 880 SUBROUTINE vlz(q,pente_max,masse,w,iq) 881 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 759 882 c 760 883 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 779 902 c Arguments: 780 903 c ---------- 781 REAL masse(ip1jmp1,llm ),pente_max782 REAL q(ip1jmp1,llm )904 REAL masse(ip1jmp1,llm,nqtot),pente_max 905 REAL q(ip1jmp1,llm,nqtot) 783 906 REAL w(ip1jmp1,llm+1) 907 INTEGER iq 784 908 c 785 909 c Local … … 792 916 REAL dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax 793 917 REAL sigw 918 919 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi 920 INTEGER ifils,iq2 ! CRisi 794 921 795 922 LOGICAL testcpu … … 805 932 c On oriente tout dans le sens de la pression c'est a dire dans le 806 933 c sens de W 934 935 write(*,*) 'vlz 923: entree' 807 936 808 937 #ifdef BIDON … … 813 942 DO l=2,llm 814 943 DO ij=1,ip1jmp1 815 dzqw(ij,l)=q(ij,l-1 )-q(ij,l)944 dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq) 816 945 adzqw(ij,l)=abs(dzqw(ij,l)) 817 946 ENDDO … … 835 964 ENDDO 836 965 966 write(*,*) 'vlz 954' 837 967 DO ij=1,ip1jmp1 838 968 dzq(ij,1)=0. … … 851 981 c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq 852 982 983 write(*,*) 'vlz 969' 853 984 DO l = 1,llm-1 854 985 do ij = 1,ip1jmp1 855 986 IF(w(ij,l+1).gt.0.) THEN 856 sigw=w(ij,l+1)/masse(ij,l+1) 857 wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1)) 987 sigw=w(ij,l+1)/masse(ij,l+1,iq) 988 wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1,iq) 989 & +0.5*(1.-sigw)*dzq(ij,l+1)) 858 990 ELSE 859 sigw=w(ij,l+1)/masse(ij,l )860 wq(ij,l+1)=w(ij,l+1)*(q(ij,l )-0.5*(1.+sigw)*dzq(ij,l))991 sigw=w(ij,l+1)/masse(ij,l,iq) 992 wq(ij,l+1)=w(ij,l+1)*(q(ij,l,iq)-0.5*(1.+sigw)*dzq(ij,l)) 861 993 ENDIF 862 994 ENDDO … … 868 1000 ENDDO 869 1001 1002 ! CRisi: appel récursif de l'advection sur les fils. 1003 ! Il faut faire ça avant d'avoir mis à jour q et masse 1004 write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq) 1005 if (nqfils(iq).gt.0) then 1006 do ifils=1,nqdesc(iq) 1007 iq2=iqfils(ifils,iq) 1008 DO l=1,llm 1009 DO ij=1,ip1jmp1 1010 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 1011 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 1012 enddo 1013 enddo 1014 enddo !do ifils=1,nqdesc(iq) 1015 1016 do ifils=1,nqfils(iq) 1017 iq2=iqfils(ifils,iq) 1018 call vlz(Ratio,pente_max,masseq,wq,iq2) 1019 enddo !do ifils=1,nqfils(iq) 1020 endif !if (nqfils(iq).gt.0) then 1021 ! end CRisi 1022 870 1023 DO l=1,llm 871 1024 DO ij=1,ip1jmp1 872 newmasse=masse(ij,l )+w(ij,l+1)-w(ij,l)873 q(ij,l )=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))1025 newmasse=masse(ij,l,iq)+w(ij,l+1)-w(ij,l) 1026 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+wq(ij,l+1)-wq(ij,l)) 874 1027 & /newmasse 875 masse(ij,l)=newmasse 876 ENDDO 877 ENDDO 878 1028 masse(ij,l,iq)=newmasse 1029 ENDDO 1030 ENDDO 1031 1032 ! retablir les fils en rapport de melange par rapport a l'air: 1033 if (nqfils(iq).gt.0) then 1034 do ifils=1,nqdesc(iq) 1035 iq2=iqfils(ifils,iq) 1036 DO l=1,llm 1037 DO ij=1,ip1jmp1 1038 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 1039 enddo 1040 enddo 1041 enddo !do ifils=1,nqdesc(iq) 1042 endif !if (nqfils(iq).gt.0) then 1043 write(*,*) 'vlsplt 1032' 879 1044 880 1045 RETURN -
LMDZ5/trunk/libf/dyn3d/vlspltqs.F
r1907 r2270 3 3 c 4 4 SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt, 5 , p,pk,teta ) 5 , p,pk,teta,iq ) 6 USE infotrac, ONLY: nqtot,nqdesc,iqfils 6 7 c 7 8 c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron … … 35 36 REAL masse(ip1jmp1,llm),pente_max 36 37 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm) 37 REAL q(ip1jmp1,llm )38 REAL q(ip1jmp1,llm,nqtot) 38 39 REAL w(ip1jmp1,llm),pdt 39 40 REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm) 41 INTEGER iq ! CRisi 40 42 c 41 43 c Local … … 43 45 c 44 46 INTEGER i,ij,l,j,ii 47 INTEGER ifils,iq2 ! CRisi 45 48 c 46 49 REAL qsat(ip1jmp1,llm) 47 REAL zm(ip1jmp1,llm )50 REAL zm(ip1jmp1,llm,nqtot) 48 51 REAL mu(ip1jmp1,llm) 49 52 REAL mv(ip1jm,llm) 50 53 REAL mw(ip1jmp1,llm+1) 51 REAL zq(ip1jmp1,llm )54 REAL zq(ip1jmp1,llm,nqtot) 52 55 REAL temps1,temps2,temps3 53 56 REAL zzpbar, zzw … … 116 119 ENDDO 117 120 118 CALL SCOPY(ijp1llm,q,1,zq,1) 119 CALL SCOPY(ijp1llm,masse,1,zm,1) 121 CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1) 122 CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1) 123 if (nqdesc(iq).gt.0) then 124 do ifils=1,nqdesc(iq) 125 iq2=iqfils(ifils,iq) 126 CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1) 127 enddo 128 endif !if (nqfils(iq).gt.0) then 120 129 121 130 c call minmaxq(zq,qmin,qmax,'avant vlxqs ') 122 call vlxqs(zq,pente_max,zm,mu,qsat) 123 131 call vlxqs(zq,pente_max,zm,mu,qsat,iq) 124 132 125 133 c call minmaxq(zq,qmin,qmax,'avant vlyqs ') 126 134 127 call vlyqs(zq,pente_max,zm,mv,qsat) 128 135 call vlyqs(zq,pente_max,zm,mv,qsat,iq) 129 136 130 137 c call minmaxq(zq,qmin,qmax,'avant vlz ') 131 138 132 call vlz(zq,pente_max,zm,mw) 133 139 call vlz(zq,pente_max,zm,mw,iq) 134 140 135 141 c call minmaxq(zq,qmin,qmax,'avant vlyqs ') 136 142 c call minmaxq(zm,qmin,qmax,'M avant vlyqs ') 137 143 138 call vlyqs(zq,pente_max,zm,mv,qsat) 139 144 call vlyqs(zq,pente_max,zm,mv,qsat,iq) 140 145 141 146 c call minmaxq(zq,qmin,qmax,'avant vlxqs ') 142 147 c call minmaxq(zm,qmin,qmax,'M avant vlxqs ') 143 148 144 call vlxqs(zq,pente_max,zm,mu,qsat )149 call vlxqs(zq,pente_max,zm,mu,qsat,iq) 145 150 146 151 c call minmaxq(zq,qmin,qmax,'apres vlxqs ') … … 150 155 DO l=1,llm 151 156 DO ij=1,ip1jmp1 152 q(ij,l )=zq(ij,l)157 q(ij,l,iq)=zq(ij,l,iq) 153 158 ENDDO 154 159 DO ij=1,ip1jm+1,iip1 155 q(ij+iim,l)=q(ij,l) 156 ENDDO 157 ENDDO 160 q(ij+iim,l,iq)=q(ij,l,iq) 161 ENDDO 162 ENDDO 163 ! CRisi: aussi pour les fils 164 if (nqdesc(iq).gt.0) then 165 do ifils=1,nqdesc(iq) 166 iq2=iqfils(ifils,iq) 167 DO l=1,llm 168 DO ij=1,ip1jmp1 169 q(ij,l,iq2)=zq(ij,l,iq2) 170 ENDDO 171 DO ij=1,ip1jm+1,iip1 172 q(ij+iim,l,iq2)=q(ij,l,iq2) 173 ENDDO 174 ENDDO 175 enddo !do ifils=1,nqdesc(iq) 176 endif ! if (nqfils(iq).gt.0) then 177 write(*,*) 'vlspltqs 183: fin de la routine' 158 178 159 179 RETURN 160 180 END 161 SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat) 181 SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq) 182 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 183 162 184 c 163 185 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 179 201 c Arguments: 180 202 c ---------- 181 REAL masse(ip1jmp1,llm ),pente_max203 REAL masse(ip1jmp1,llm,nqtot),pente_max 182 204 REAL u_m( ip1jmp1,llm ) 183 REAL q(ip1jmp1,llm )205 REAL q(ip1jmp1,llm,nqtot) 184 206 REAL qsat(ip1jmp1,llm) 207 INTEGER iq ! CRisi 185 208 c 186 209 c Local … … 195 218 REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm) 196 219 REAL u_mq(ip1jmp1,llm) 220 221 ! CRisi 222 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) 223 INTEGER ifils,iq2 ! CRisi 197 224 198 225 Logical first,testcpu … … 227 254 DO l = 1, llm 228 255 DO ij=iip2,ip1jm-1 229 dxqu(ij)=q(ij+1,l )-q(ij,l)256 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 230 257 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0' 231 c sigu(ij)=u_m(ij,l)/masse(ij,l )258 c sigu(ij)=u_m(ij,l)/masse(ij,l,iq) 232 259 ENDDO 233 260 DO ij=iip1+iip1,ip1jm,iip1 … … 281 308 DO l = 1, llm 282 309 DO ij=iip2,ip1jm-1 283 dxqu(ij)=q(ij+1,l )-q(ij,l)310 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 284 311 ENDDO 285 312 DO ij=iip1+iip1,ip1jm,iip1 … … 323 350 DO l=1,llm 324 351 DO ij=iip2,ip1jm-1 325 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l ),326 , 1.+u_m(ij,l)/masse(ij+1,l ),352 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq), 353 , 1.+u_m(ij,l)/masse(ij+1,l,iq), 327 354 , u_m(ij,l)) 328 355 zdum(ij,l)=0.5*zdum(ij,l) 329 356 u_mq(ij,l)=cvmgp( 330 , q(ij,l )+zdum(ij,l)*dxq(ij,l),331 , q(ij+1,l )-zdum(ij,l)*dxq(ij+1,l),357 , q(ij,l,iq)+zdum(ij,l)*dxq(ij,l), 358 , q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l), 332 359 , u_m(ij,l)) 333 360 u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l) … … 341 368 DO ij=iip2,ip1jm-1 342 369 IF (u_m(ij,l).gt.0.) THEN 343 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l )370 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq) 344 371 u_mq(ij,l)=u_m(ij,l)* 345 $ min(q(ij,l )+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))372 $ min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l)) 346 373 ELSE 347 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l )374 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq) 348 375 u_mq(ij,l)=u_m(ij,l)* 349 $ min(q(ij+1,l )-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))376 $ min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l)) 350 377 ENDIF 351 378 ENDDO … … 416 443 i=ijq-(j-1)*iip1 417 444 c accumulation pour les mailles completements advectees 418 do while(zu_m.gt.masse(ijq,l)) 419 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l) 420 zu_m=zu_m-masse(ijq,l) 445 do while(zu_m.gt.masse(ijq,l,iq)) 446 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq) 447 & *masse(ijq,l,iq) 448 zu_m=zu_m-masse(ijq,l,iq) 421 449 i=mod(i-2+iim,iim)+1 422 450 ijq=(j-1)*iip1+i … … 424 452 c ajout de la maille non completement advectee 425 453 u_mq(ij,l)=u_mq(ij,l)+zu_m* 426 & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l)) 454 & (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq)) 455 & *dxq(ijq,l)) 427 456 ELSE 428 457 ijq=ij+1 429 458 i=ijq-(j-1)*iip1 430 459 c accumulation pour les mailles completements advectees 431 do while(-zu_m.gt.masse(ijq,l)) 432 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l) 433 zu_m=zu_m+masse(ijq,l) 460 do while(-zu_m.gt.masse(ijq,l,iq)) 461 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq) 462 & *masse(ijq,l,iq) 463 zu_m=zu_m+masse(ijq,l,iq) 434 464 i=mod(i,iim)+1 435 465 ijq=(j-1)*iip1+i 436 466 ENDDO 437 467 c ajout de la maille non completement advectee 438 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l )-439 & 0.5*(1.+zu_m/masse(ijq,l ))*dxq(ijq,l))468 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)- 469 & 0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 440 470 ENDIF 441 471 ENDDO … … 454 484 ENDDO 455 485 486 ! CRisi: appel récursif de l'advection sur les fils. 487 ! Il faut faire ça avant d'avoir mis à jour q et masse 488 write(*,*) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq) 489 490 if (nqfils(iq).gt.0) then 491 do ifils=1,nqdesc(iq) 492 iq2=iqfils(ifils,iq) 493 DO l=1,llm 494 DO ij=iip2,ip1jm 495 ! On a besoin de q et masse seulement entre iip2 et ip1jm 496 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 497 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 498 enddo 499 enddo 500 enddo !do ifils=1,nqdesc(iq) 501 do ifils=1,nqfils(iq) 502 iq2=iqfils(ifils,iq) 503 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 504 enddo !do ifils=1,nqfils(iq) 505 endif !if (nqfils(iq).gt.0) then 506 ! end CRisi 456 507 457 508 c calcul des tendances … … 459 510 DO l=1,llm 460 511 DO ij=iip2+1,ip1jm 461 new_m=masse(ij,l )+u_m(ij-1,l)-u_m(ij,l)462 q(ij,l )=(q(ij,l)*masse(ij,l)+512 new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l) 513 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+ 463 514 & u_mq(ij-1,l)-u_mq(ij,l)) 464 515 & /new_m 465 masse(ij,l )=new_m516 masse(ij,l,iq)=new_m 466 517 ENDDO 467 518 c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous) 468 519 DO ij=iip1+iip1,ip1jm,iip1 469 q(ij-iim,l)=q(ij,l) 470 masse(ij-iim,l)=masse(ij,l) 471 ENDDO 472 ENDDO 520 q(ij-iim,l,iq)=q(ij,l,iq) 521 masse(ij-iim,l,iq)=masse(ij,l,iq) 522 ENDDO 523 ENDDO 524 525 ! retablir les fils en rapport de melange par rapport a l'air: 526 ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio 527 ! puis on boucle en longitude 528 if (nqdesc(iq).gt.0) then 529 do ifils=1,nqdesc(iq) 530 iq2=iqfils(ifils,iq) 531 DO l=1,llm 532 DO ij=iip2+1,ip1jm 533 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 534 enddo 535 DO ij=iip1+iip1,ip1jm,iip1 536 q(ij-iim,l,iq2)=q(ij,l,iq2) 537 enddo ! DO ij=ijb+iip1-1,ije,iip1 538 enddo !DO l=1,llm 539 enddo !do ifils=1,nqdesc(iq) 540 endif !if (nqfils(iq).gt.0) then 473 541 474 542 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) … … 478 546 RETURN 479 547 END 480 SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat) 548 SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq) 549 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 481 550 c 482 551 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 502 571 c Arguments: 503 572 c ---------- 504 REAL masse(ip1jmp1,llm ),pente_max573 REAL masse(ip1jmp1,llm,nqtot),pente_max 505 574 REAL masse_adv_v( ip1jm,llm) 506 REAL q(ip1jmp1,llm )575 REAL q(ip1jmp1,llm,nqtot) 507 576 REAL qsat(ip1jmp1,llm) 577 INTEGER iq ! CRisi 508 578 c 509 579 c Local … … 529 599 SAVE sinlon,coslon,sinlondlon,coslondlon 530 600 SAVE airej2,airejjm 601 602 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi 603 INTEGER ifils,iq2 ! CRisi 531 604 c 532 605 c … … 567 640 568 641 DO i = 1, iim 569 airescb(i) = aire(i+ iip1) * q(i+ iip1,l )570 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l )642 airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq) 643 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq) 571 644 ENDDO 572 645 qpns = SSUM( iim, airescb ,1 ) / airej2 … … 576 649 577 650 DO ij=1,ip1jm 578 dyqv(ij)=q(ij,l )-q(ij+iip1,l)651 dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq) 579 652 adyqv(ij)=abs(dyqv(ij)) 580 653 ENDDO … … 591 664 592 665 DO ij=1,iip1 593 dyq(ij,l)=qpns-q(ij+iip1,l )594 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l )-qpsn666 dyq(ij,l)=qpns-q(ij+iip1,l,iq) 667 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn 595 668 ENDDO 596 669 … … 710 783 DO ij=1,ip1jm 711 784 IF( masse_adv_v(ij,l).GT.0. ) THEN 712 qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l ) + 713 , dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))) 785 qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq ) + 786 , dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l) 787 , /masse(ij+iip1,l,iq))) 714 788 ELSE 715 qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l ) - dyq(ij,l) *716 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l )) )789 qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) * 790 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) ) 717 791 ENDIF 718 792 qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l) … … 721 795 722 796 797 ! CRisi: appel récursif de l'advection sur les fils. 798 ! Il faut faire ça avant d'avoir mis à jour q et masse 799 write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq) 800 801 if (nqfils(iq).gt.0) then 802 do ifils=1,nqdesc(iq) 803 iq2=iqfils(ifils,iq) 804 DO l=1,llm 805 DO ij=1,ip1jmp1 806 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 807 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 808 enddo 809 enddo 810 enddo !do ifils=1,nqdesc(iq) 811 812 do ifils=1,nqfils(iq) 813 iq2=iqfils(ifils,iq) 814 write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2 815 call vly(Ratio,pente_max,masseq,qbyv,iq2) 816 enddo !do ifils=1,nqfils(iq) 817 endif !if (nqfils(iq).gt.0) then 818 723 819 DO l=1,llm 724 820 DO ij=iip2,ip1jm 725 newmasse=masse(ij,l )821 newmasse=masse(ij,l,iq) 726 822 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l) 727 q(ij,l )=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))728 & /newmasse729 masse(ij,l )=newmasse823 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l) 824 & -qbyv(ij-iip1,l))/newmasse 825 masse(ij,l,iq)=newmasse 730 826 ENDDO 731 827 c.-. ancienne version … … 733 829 convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln 734 830 DO ij = 1,iip1 735 newmasse=masse(ij,l )+convmpn*aire(ij)736 q(ij,l )=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/831 newmasse=masse(ij,l,iq)+convmpn*aire(ij) 832 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/ 737 833 & newmasse 738 masse(ij,l )=newmasse834 masse(ij,l,iq)=newmasse 739 835 ENDDO 740 836 convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols 741 837 convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols 742 838 DO ij = ip1jm+1,ip1jmp1 743 newmasse=masse(ij,l )+convmps*aire(ij)744 q(ij,l )=(q(ij,l)*masse(ij,l)+convps*aire(ij))/839 newmasse=masse(ij,l,iq)+convmps*aire(ij) 840 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/ 745 841 & newmasse 746 masse(ij,l )=newmasse842 masse(ij,l,iq)=newmasse 747 843 ENDDO 748 844 c.-. fin ancienne version … … 757 853 c DO ij = 1,iip1 758 854 c q(ij,l)=newq 759 c masse(ij,l )=newmasse*aire(ij)855 c masse(ij,l,iq)=newmasse*aire(ij) 760 856 c ENDDO 761 857 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) … … 767 863 c DO ij = ip1jm+1,ip1jmp1 768 864 c q(ij,l)=newq 769 c masse(ij,l )=newmasse*aire(ij)865 c masse(ij,l,iq)=newmasse*aire(ij) 770 866 c ENDDO 771 867 c._. fin nouvelle version 772 868 ENDDO 773 869 870 write(*,*) 'vly 866' 871 872 ! retablir les fils en rapport de melange par rapport a l'air: 873 if (nqdesc(iq).gt.0) then 874 do ifils=1,nqdesc(iq) 875 iq2=iqfils(ifils,iq) 876 DO l=1,llm 877 DO ij=1,ip1jmp1 878 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 879 enddo 880 enddo 881 enddo !do ifils=1,nqdesc(iq) 882 endif !if (nqfils(iq).gt.0) then 883 write(*,*) 'vly 879' 884 774 885 RETURN 775 886 END -
LMDZ5/trunk/libf/dyn3d_common/infotrac.F90
r2262 r2270 12 12 INTEGER, SAVE :: nbtr 13 13 14 ! CRisi: nb traceurs pères= directement advectés par l'air 15 INTEGER, SAVE :: nqperes 16 14 17 ! Name variables 15 18 CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics … … 22 25 ! dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code. 23 26 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: niadv ! equivalent dyn / physique 27 28 ! CRisi: tableaux de fils 29 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqfils 30 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqdesc ! nombres de fils + nombre de tous les petits fils sur toutes les générations 31 INTEGER, SAVE :: nqdesc_tot 32 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqfils 33 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iqpere 24 34 25 35 ! conv_flg(it)=0 : convection desactivated for tracer number it … … 30 40 CHARACTER(len=4),SAVE :: type_trac 31 41 CHARACTER(len=8),DIMENSION(:),ALLOCATABLE, SAVE :: solsym 42 43 ! CRisi: cas particulier des isotopes 44 LOGICAL,SAVE :: ok_isotopes,ok_iso_verif,ok_isotrac,ok_init_iso 45 INTEGER :: niso_possibles 46 PARAMETER ( niso_possibles=5) 47 real, DIMENSION (niso_possibles),SAVE :: tnat,alpha_ideal 48 LOGICAL, DIMENSION(niso_possibles),SAVE :: use_iso 49 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqiso ! donne indice iq en fn de (ixt,phase) 50 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_num ! donne numéro iso entre 1 et niso_possibles en fn de nqtot 51 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_indnum ! donne numéro iso entre 1 et niso effectif en fn de nqtot 52 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: zone_num ! donne numéro de la zone de tracage en fn de nqtot 53 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: phase_num ! donne numéro de la zone de tracage en fn de nqtot 54 INTEGER, DIMENSION(niso_possibles), SAVE :: indnum_fn_num ! donne indice entre entre 1 et niso en fonction du numéro d isotope entre 1 et niso_possibles 55 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: index_trac ! numéro ixt en fn izone, indnum entre 1 et niso 56 INTEGER,SAVE :: niso,ntraceurs_zone,ntraciso 32 57 33 58 CONTAINS … … 63 88 64 89 CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_0 ! tracer short name 90 CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name: CRisi 65 91 CHARACTER(len=3), DIMENSION(30) :: descrq 66 92 CHARACTER(len=1), DIMENSION(3) :: txts … … 70 96 INTEGER :: nqtrue ! number of tracers read from tracer.def, without higer order of moment 71 97 INTEGER :: iq, new_iq, iiq, jq, ierr 98 INTEGER :: ifils,ipere,generation ! CRisi 99 LOGICAL :: continu,nouveau_traceurdef 100 INTEGER :: IOstatus ! gestion de la retrocompatibilite de traceur.def 101 CHARACTER(len=15) :: tchaine 72 102 73 103 character(len=*),parameter :: modname="infotrac_init" … … 134 164 WRITE(lunout,*) trim(modname),': Open traceur.def : ok' 135 165 READ(90,*) nqtrue 166 write(lunout,*) 'nqtrue=',nqtrue 136 167 ELSE 137 168 WRITE(lunout,*) trim(modname),': Problem in opening traceur.def' … … 192 223 ! Allocate variables depending on nqtrue 193 224 ! 194 ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue) )225 ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue),tnom_transp(nqtrue)) 195 226 ! 196 227 !jyg< … … 230 261 ! Continue to read tracer.def 231 262 DO iq=1,nqtrue 232 READ(90,*) hadv(iq),vadv(iq),tnom_0(iq) 233 END DO 263 264 write(*,*) 'infotrac 237: iq=',iq 265 ! CRisi: ajout du nom du fluide transporteur 266 ! mais rester retro compatible 267 READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine 268 write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq) 269 write(lunout,*) 'tchaine=',trim(tchaine) 270 write(*,*) 'infotrac 238: IOstatus=',IOstatus 271 if (IOstatus.ne.0) then 272 CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1) 273 endif 274 ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un 275 ! espace ou pas au milieu de la chaine. 276 continu=1 277 nouveau_traceurdef=0 278 iiq=1 279 do while (continu) 280 if (tchaine(iiq:iiq).eq.' ') then 281 nouveau_traceurdef=1 282 continu=0 283 else if (iiq.lt.LEN_TRIM(tchaine)) then 284 iiq=iiq+1 285 else 286 continu=0 287 endif 288 enddo 289 write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef 290 if (nouveau_traceurdef) then 291 write(lunout,*) 'C''est la nouvelle version de traceur.def' 292 tnom_0(iq)=tchaine(1:iiq-1) 293 tnom_transp(iq)=tchaine(iiq+1:15) 294 else 295 write(lunout,*) 'C''est l''ancienne version de traceur.def' 296 write(lunout,*) 'On suppose que les traceurs sont tous d''air' 297 tnom_0(iq)=tchaine 298 tnom_transp(iq) = 'air' 299 endif 300 write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>' 301 write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>' 302 303 END DO !DO iq=1,nqtrue 234 304 CLOSE(90) 305 235 306 ELSE ! Without tracer.def, set default values 236 307 if (planet_type=="earth") then … … 239 310 vadv(1) = 14 240 311 tnom_0(1) = 'H2Ov' 312 tnom_transp(1) = 'air' 241 313 hadv(2) = 10 242 314 vadv(2) = 10 243 315 tnom_0(2) = 'H2Ol' 316 tnom_transp(2) = 'air' 244 317 hadv(3) = 10 245 318 vadv(3) = 10 246 319 tnom_0(3) = 'RN' 320 tnom_transp(3) = 'air' 247 321 hadv(4) = 10 248 322 vadv(4) = 10 249 323 tnom_0(4) = 'PB' 324 tnom_transp(4) = 'air' 250 325 else ! default for other planets 251 326 hadv(1) = 10 252 327 vadv(1) = 10 253 328 tnom_0(1) = 'dummy' 329 tnom_transp(1) = 'dummy' 254 330 endif ! of if (planet_type=="earth") 255 331 END IF … … 258 334 WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue 259 335 DO iq=1,nqtrue 260 WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq) 336 WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq),tnom_transp(iq) 261 337 END DO 262 338 … … 447 523 END DO 448 524 525 526 ! CRisi: quels sont les traceurs fils et les traceurs pères. 527 ! initialiser tous les tableaux d'indices liés aux traceurs familiaux 528 ! + vérifier que tous les pères sont écrits en premières positions 529 ALLOCATE(nqfils(nqtot),nqdesc(nqtot)) 530 ALLOCATE(iqfils(nqtot,nqtot)) 531 ALLOCATE(iqpere(nqtot)) 532 nqperes=0 533 nqfils(:)=0 534 nqdesc(:)=0 535 iqfils(:,:)=0 536 iqpere(:)=0 537 nqdesc_tot=0 538 DO iq=1,nqtot 539 if (tnom_transp(iq) == 'air') then 540 ! ceci est un traceur père 541 WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un pere' 542 nqperes=nqperes+1 543 iqpere(iq)=0 544 else !if (tnom_transp(iq) == 'air') then 545 ! ceci est un fils. Qui est son père? 546 WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un fils' 547 continu=.true. 548 ipere=1 549 do while (continu) 550 if (tnom_transp(iq) == tnom_0(ipere)) then 551 ! Son père est ipere 552 WRITE(lunout,*) 'Le traceur',iq,'appele ', & 553 & trim(tnom_0(iq)),' est le fils de ',ipere,'appele ',trim(tnom_0(ipere)) 554 nqfils(ipere)=nqfils(ipere)+1 555 iqfils(nqfils(ipere),ipere)=iq 556 iqpere(iq)=ipere 557 continu=.false. 558 else !if (tnom_transp(iq) == tnom_0(ipere)) then 559 ipere=ipere+1 560 if (ipere.gt.nqtot) then 561 WRITE(lunout,*) 'Le traceur',iq,'appele ', & 562 & trim(tnom_0(iq)),', est orpelin.' 563 CALL abort_gcm('infotrac_init','Un traceur est orphelin',1) 564 endif !if (ipere.gt.nqtot) then 565 endif !if (tnom_transp(iq) == tnom_0(ipere)) then 566 enddo !do while (continu) 567 endif !if (tnom_transp(iq) == 'air') then 568 enddo !DO iq=1,nqtot 569 WRITE(lunout,*) 'infotrac: nqperes=',nqperes 570 WRITE(lunout,*) 'nqfils=',nqfils 571 WRITE(lunout,*) 'iqpere=',iqpere 572 WRITE(lunout,*) 'iqfils=',iqfils 573 574 ! Calculer le nombre de descendants à partir de iqfils et de nbfils 575 DO iq=1,nqtot 576 generation=0 577 continu=.true. 578 ifils=iq 579 do while (continu) 580 ipere=iqpere(ifils) 581 if (ipere.gt.0) then 582 nqdesc(ipere)=nqdesc(ipere)+1 583 nqdesc_tot=nqdesc_tot+1 584 iqfils(nqdesc(ipere),ipere)=iq 585 ifils=ipere 586 generation=generation+1 587 else !if (ipere.gt.0) then 588 continu=.false. 589 endif !if (ipere.gt.0) then 590 enddo !do while (continu) 591 WRITE(lunout,*) 'Le traceur ',iq,', appele ',trim(tnom_0(iq)),' est un traceur de generation: ',generation 592 enddo !DO iq=1,nqtot 593 WRITE(lunout,*) 'infotrac: nqdesc=',nqdesc 594 WRITE(lunout,*) 'iqfils=',iqfils 595 WRITE(lunout,*) 'nqdesc_tot=',nqdesc_tot 596 597 ! Interdire autres schémas que 10 pour les traceurs fils, et autres schémas 598 ! que 10 et 14 si des pères ont des fils 599 do iq=1,nqtot 600 if (iqpere(iq).gt.0) then 601 ! ce traceur a un père qui n'est pas l'air 602 ! Seul le schéma 10 est autorisé 603 if (iadv(iq)/=10) then 604 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for sons' 605 CALL abort_gcm('infotrac_init','Sons should be advected by scheme 10',1) 606 endif 607 ! Le traceur père ne peut être advecté que par schéma 10 ou 14: 608 IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN 609 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for fathers' 610 CALL abort_gcm('infotrac_init','Fathers should be advected by scheme 10 ou 14',1) 611 endif !IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN 612 endif !if (iqpere(iq).gt.0) the 613 enddo !do iq=1,nqtot 614 615 616 ! detecter quels sont les traceurs isotopiques parmi des traceurs 617 call infotrac_isoinit(tnom_0,nqtrue) 618 449 619 !----------------------------------------------------------------------- 450 620 ! Finalize : 451 621 ! 452 DEALLOCATE(tnom_0, hadv, vadv )622 DEALLOCATE(tnom_0, hadv, vadv,tnom_transp) 453 623 454 624 455 625 END SUBROUTINE infotrac_init 456 626 627 SUBROUTINE infotrac_isoinit(tnom_0,nqtrue) 628 629 #ifdef CPP_IOIPSL 630 use IOIPSL 631 #else 632 ! if not using IOIPSL, we still need to use (a local version of) getin 633 use ioipsl_getincom 634 #endif 635 implicit none 636 637 ! inputs 638 INTEGER nqtrue 639 CHARACTER(len=15) tnom_0(nqtrue) 640 641 ! locals 642 CHARACTER(len=3), DIMENSION(niso_possibles) :: tnom_iso 643 INTEGER, ALLOCATABLE,DIMENSION(:,:) :: nb_iso,nb_traciso 644 INTEGER, ALLOCATABLE,DIMENSION(:) :: nb_isoind 645 INTEGER :: ntraceurs_zone_prec,iq,phase,ixt,iiso,izone 646 CHARACTER(len=19) :: tnom_trac 647 INCLUDE "iniprint.h" 648 649 tnom_iso=(/'eau','HDO','O18','O17','HTO'/) 650 651 ALLOCATE(nb_iso(niso_possibles,nqo)) 652 ALLOCATE(nb_isoind(nqo)) 653 ALLOCATE(nb_traciso(niso_possibles,nqo)) 654 ALLOCATE(iso_num(nqtot)) 655 ALLOCATE(iso_indnum(nqtot)) 656 ALLOCATE(zone_num(nqtot)) 657 ALLOCATE(phase_num(nqtot)) 658 659 iso_num(:)=0 660 iso_indnum(:)=0 661 zone_num(:)=0 662 phase_num(:)=0 663 indnum_fn_num(:)=0 664 use_iso(:)=.false. 665 nb_iso(:,:)=0 666 nb_isoind(:)=0 667 nb_traciso(:,:)=0 668 niso=0 669 ntraceurs_zone=0 670 ntraceurs_zone_prec=0 671 ntraciso=0 672 673 do iq=nqo+1,nqtot 674 write(lunout,*) 'infotrac 569: iq,tnom_0(iq)=',iq,tnom_0(iq) 675 do phase=1,nqo 676 do ixt= 1,niso_possibles 677 tnom_trac=trim(tnom_0(phase))//'_' 678 tnom_trac=trim(tnom_trac)//trim(tnom_iso(ixt)) 679 write(*,*) 'phase,ixt,tnom_trac=',phase,ixt,tnom_trac 680 IF (tnom_0(iq) == tnom_trac) then 681 write(lunout,*) 'Ce traceur est un isotope' 682 nb_iso(ixt,phase)=nb_iso(ixt,phase)+1 683 nb_isoind(phase)=nb_isoind(phase)+1 684 iso_num(iq)=ixt 685 iso_indnum(iq)=nb_isoind(phase) 686 indnum_fn_num(ixt)=iso_indnum(iq) 687 phase_num(iq)=phase 688 write(lunout,*) 'iso_num(iq)=',iso_num(iq) 689 write(lunout,*) 'iso_indnum(iq)=',iso_indnum(iq) 690 write(lunout,*) 'indnum_fn_num(ixt)=',indnum_fn_num(ixt) 691 write(lunout,*) 'phase_num(iq)=',phase_num(iq) 692 goto 20 693 else if (iqpere(iq).gt.0) then 694 if (tnom_0(iqpere(iq)) == tnom_trac) then 695 write(lunout,*) 'Ce traceur est le fils d''un isotope' 696 ! c'est un traceur d'isotope 697 nb_traciso(ixt,phase)=nb_traciso(ixt,phase)+1 698 iso_num(iq)=ixt 699 iso_indnum(iq)=indnum_fn_num(ixt) 700 zone_num(iq)=nb_traciso(ixt,phase) 701 phase_num(iq)=phase 702 write(lunout,*) 'iso_num(iq)=',iso_num(iq) 703 write(lunout,*) 'phase_num(iq)=',phase_num(iq) 704 write(lunout,*) 'zone_num(iq)=',zone_num(iq) 705 goto 20 706 endif !if (tnom_0(iqpere(iq)) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then 707 endif !IF (tnom_0(iq) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then 708 enddo !do ixt= niso_possibles 709 enddo !do phase=1,nqo 710 20 continue 711 enddo !do iq=1,nqtot 712 713 write(lunout,*) 'iso_num=',iso_num 714 write(lunout,*) 'iso_indnum=',iso_indnum 715 write(lunout,*) 'zone_num=',zone_num 716 write(lunout,*) 'phase_num=',phase_num 717 write(lunout,*) 'indnum_fn_num=',indnum_fn_num 718 719 do ixt= 1,niso_possibles 720 721 if (nb_iso(ixt,1).eq.1) then 722 ! on vérifie que toutes les phases ont le même nombre de 723 ! traceurs 724 do phase=2,nqo 725 if (nb_iso(ixt,phase).ne.nb_iso(ixt,1)) then 726 write(lunout,*) 'ixt,phase,nb_iso=',ixt,phase,nb_iso(ixt,phase) 727 CALL abort_gcm('infotrac_init','Phases must have same number of isotopes',1) 728 endif 729 enddo !do phase=2,nqo 730 731 niso=niso+1 732 use_iso(ixt)=.true. 733 ntraceurs_zone=nb_traciso(ixt,1) 734 735 ! on vérifie que toutes les phases ont le même nombre de 736 ! traceurs 737 do phase=2,nqo 738 if (nb_traciso(ixt,phase).ne.ntraceurs_zone) then 739 write(lunout,*) 'ixt,phase,nb_traciso=',ixt,phase,nb_traciso(ixt,phase) 740 write(lunout,*) 'ntraceurs_zone=',ntraceurs_zone 741 CALL abort_gcm('infotrac_init','Phases must have same number of tracers',1) 742 endif 743 enddo !do phase=2,nqo 744 ! on vérifie que tous les isotopes ont le même nombre de 745 ! traceurs 746 if (ntraceurs_zone_prec.gt.0) then 747 if (ntraceurs_zone.eq.ntraceurs_zone_prec) then 748 ntraceurs_zone_prec=ntraceurs_zone 749 else !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then 750 write(*,*) 'ntraceurs_zone_prec,ntraceurs_zone=',ntraceurs_zone_prec,ntraceurs_zone 751 CALL abort_gcm('infotrac_init', & 752 &'Isotope tracers are not well defined in traceur.def',1) 753 endif !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then 754 endif !if (ntraceurs_zone_prec.gt.0) then 755 756 else if (nb_iso(ixt,1).ne.0) then 757 WRITE(lunout,*) 'nqo,ixt=',nqo,ixt 758 WRITE(lunout,*) 'nb_iso(ixt,1)=',nb_iso(ixt,1) 759 CALL abort_gcm('infotrac_init','Isotopes are not well defined in traceur.def',1) 760 endif !if (nb_iso(ixt,1).eq.1) then 761 enddo ! do ixt= niso_possibles 762 763 ! dimensions isotopique: 764 ntraciso=niso*(ntraceurs_zone+1) 765 WRITE(lunout,*) 'niso=',niso 766 WRITE(lunout,*) 'ntraceurs_zone,ntraciso=',ntraceurs_zone,ntraciso 767 768 ! flags isotopiques: 769 if (niso.gt.0) then 770 ok_isotopes=.true. 771 else 772 ok_isotopes=.false. 773 endif 774 WRITE(lunout,*) 'ok_isotopes=',ok_isotopes 775 776 if (ok_isotopes) then 777 ok_iso_verif=.false. 778 call getin('ok_iso_verif',ok_iso_verif) 779 ok_init_iso=.false. 780 call getin('ok_init_iso',ok_init_iso) 781 tnat=(/1.0,155.76e-6,2005.2e-6,0.004/100.,0.0/) 782 alpha_ideal=(/1.0,1.01,1.006,1.003,1.0/) 783 endif !if (ok_isotopes) then 784 WRITE(lunout,*) 'ok_iso_verif=',ok_iso_verif 785 WRITE(lunout,*) 'ok_init_iso=',ok_init_iso 786 787 if (ntraceurs_zone.gt.0) then 788 ok_isotrac=.true. 789 else 790 ok_isotrac=.false. 791 endif 792 WRITE(lunout,*) 'ok_isotrac=',ok_isotrac 793 794 ! remplissage du tableau iqiso(ntraciso,phase) 795 ALLOCATE(iqiso(ntraciso,nqo)) 796 iqiso(:,:)=0 797 do iq=1,nqtot 798 if (iso_num(iq).gt.0) then 799 ixt=iso_indnum(iq)+zone_num(iq)*niso 800 iqiso(ixt,phase_num(iq))=iq 801 endif 802 enddo 803 WRITE(lunout,*) 'iqiso=',iqiso 804 805 ! replissage du tableau index_trac(ntraceurs_zone,niso) 806 ALLOCATE(index_trac(ntraceurs_zone,niso)) 807 if (ok_isotrac) then 808 do iiso=1,niso 809 do izone=1,ntraceurs_zone 810 index_trac(izone,iiso)=iiso+izone*niso 811 enddo 812 enddo 813 else !if (ok_isotrac) then 814 index_trac(:,:)=0.0 815 endif !if (ok_isotrac) then 816 write(lunout,*) 'index_trac=',index_trac 817 818 ! Finalize : 819 DEALLOCATE(nb_iso) 820 821 END SUBROUTINE infotrac_isoinit 822 457 823 END MODULE infotrac -
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.
