Changeset 4050 for LMDZ6/trunk/libf
- Timestamp:
- Dec 23, 2021, 6:54:17 PM (3 years ago)
- Location:
- LMDZ6/trunk/libf
- Files:
-
- 38 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/dyn3d/advtrac.F90
r2622 r4050 9 9 ! M.A Filiberti (04/2002) 10 10 ! 11 USE infotrac, ONLY: nqtot, iadv,nqperes,ok_iso_verif11 USE infotrac, ONLY: nqtot, tracers, nqperes,ok_iso_verif 12 12 USE control_mod, ONLY: iapp_tracvl, day_step 13 13 USE comconst_mod, ONLY: dtvr … … 72 72 real cflz(ip1jmp1,llm) 73 73 real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm) 74 INTEGER :: iadv 74 75 75 76 IF(iadvtr.EQ.0) THEN … … 226 227 do iq=1,nqperes 227 228 ! call clock(t_initial) 228 if(iadv(iq) == 0) cycle 229 iadv = tracers(iq)%iadv 230 SELECT CASE(iadv) 231 CASE(0); CYCLE 232 CASE(10) 229 233 ! ---------------------------------------------------------------- 230 234 ! Schema de Van Leer I MUSCL 231 235 ! ---------------------------------------------------------------- 232 if(iadv(iq).eq.10) THEN233 236 ! CRisi: on fait passer tout q pour avoir acces aux fils 234 237 … … 236 239 call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq) 237 240 241 CASE(14) 238 242 ! ---------------------------------------------------------------- 239 243 ! Schema "pseudo amont" + test sur humidite specifique 240 244 ! pour la vapeur d'eau. F. Codron 241 245 ! ---------------------------------------------------------------- 242 else if(iadv(iq).eq.14) then243 246 ! 244 247 !write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:) … … 246 249 pbarug,pbarvg,dtvr,p,pk,teta,iq) 247 250 251 CASE(12) 248 252 ! ---------------------------------------------------------------- 249 253 ! Schema de Frederic Hourdin 250 254 ! ---------------------------------------------------------------- 251 else if(iadv(iq).eq.12) then252 255 ! Pas de temps adaptatif 253 call adaptdt(iadv (iq),dtbon,n,pbarug,massem)256 call adaptdt(iadv,dtbon,n,pbarug,massem) 254 257 if (n.GT.1) then 255 258 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & … … 259 262 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1) 260 263 end do 261 else if(iadv(iq).eq.13) then264 CASE(13) 262 265 ! Pas de temps adaptatif 263 call adaptdt(iadv (iq),dtbon,n,pbarug,massem)266 call adaptdt(iadv,dtbon,n,pbarug,massem) 264 267 if (n.GT.1) then 265 268 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & … … 269 272 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2) 270 273 end do 274 CASE(20) 271 275 ! ---------------------------------------------------------------- 272 276 ! Schema de pente SLOPES 273 277 ! ---------------------------------------------------------------- 274 else if (iadv(iq).eq.20) then275 278 call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0) 276 279 280 CASE(30) 277 281 ! ---------------------------------------------------------------- 278 282 ! Schema de Prather 279 283 ! ---------------------------------------------------------------- 280 else if (iadv(iq).eq.30) then281 284 ! Pas de temps adaptatif 282 call adaptdt(iadv (iq),dtbon,n,pbarug,massem)285 call adaptdt(iadv,dtbon,n,pbarug,massem) 283 286 if (n.GT.1) then 284 287 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & … … 288 291 n,dtbon) 289 292 293 CASE(11,16,17,18) 290 294 ! ---------------------------------------------------------------- 291 295 ! Schemas PPM Lin et Rood 292 296 ! ---------------------------------------------------------------- 293 else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &294 iadv(iq).LE.18)) then295 297 296 298 ! Test sur le flux horizontal 297 299 ! Pas de temps adaptatif 298 call adaptdt(iadv (iq),dtbon,n,pbarug,massem)300 call adaptdt(iadv,dtbon,n,pbarug,massem) 299 301 if (n.GT.1) then 300 302 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & … … 316 318 317 319 !----------------------------------------------------------- 318 ! Ss-prg interface LMDZ.4->PPM3d 320 ! Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin) 319 321 !----------------------------------------------------------- 320 322 … … 327 329 ! VL (version PPM) horiz. et PPM vert. 328 330 !---------------------------------------------------------------- 329 if (iadv(iq).eq.11) then330 ! Ss-prg PPM3d de Lin331 SELECT CASE(iadv) 332 CASE(11) 331 333 call ppm3d(1,qppm(1,1,iq), & 332 334 psppm,psppm, & … … 335 337 fill,dum,220.) 336 338 339 CASE(16) 337 340 !------------------------------------------------------------- 338 341 ! Monotonic PPM 339 342 !------------------------------------------------------------- 340 else if (iadv(iq).eq.16) then341 ! Ss-prg PPM3d de Lin342 343 call ppm3d(1,qppm(1,1,iq), & 343 344 psppm,psppm, & … … 347 348 !------------------------------------------------------------- 348 349 350 CASE(17) 349 351 !------------------------------------------------------------- 350 352 ! Semi Monotonic PPM 351 353 !------------------------------------------------------------- 352 else if (iadv(iq).eq.17) then353 ! Ss-prg PPM3d de Lin354 354 call ppm3d(1,qppm(1,1,iq), & 355 355 psppm,psppm, & … … 359 359 !------------------------------------------------------------- 360 360 361 CASE(18) 361 362 !------------------------------------------------------------- 362 363 ! Positive Definite PPM 363 364 !------------------------------------------------------------- 364 else if (iadv(iq).eq.18) then365 ! Ss-prg PPM3d de Lin366 365 call ppm3d(1,qppm(1,1,iq), & 367 366 psppm,psppm, & … … 370 369 fill,dum,220.) 371 370 !------------------------------------------------------------- 372 endif371 END SELECT 373 372 enddo 374 373 !----------------------------------------------------------------- … … 376 375 !----------------------------------------------------------------- 377 376 call interpost(q(1,1,iq),qppm(1,1,iq)) 378 endif377 END SELECT 379 378 !---------------------------------------------------------------------- 380 379 -
LMDZ6/trunk/libf/dyn3d/check_isotopes.F
r4037 r4050 1 1 subroutine check_isotopes_seq(q,ip1jmp1,err_msg) 2 USE infotrac 2 USE infotrac, ONLY: nqtot, nqo, niso, ntraciso, ntraceurs_zone, 3 & ok_isotopes, ok_isotrac, use_iso, 4 & iqiso, index_trac,indnum_fn_num, tnat 3 5 implicit none 4 6 -
LMDZ6/trunk/libf/dyn3d/dynetat0.f90
r4046 r4050 6 6 ! Purpose: Initial state reading. 7 7 !------------------------------------------------------------------------------- 8 USE infotrac 8 USE infotrac, ONLY: nqtot, tracers, iqiso, iso_indnum, tnat, alpha_ideal, & 9 ok_isotopes, maxlen 9 10 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_INQ_VARID, NF90_NoErr, & 10 11 NF90_CLOSE, NF90_GET_VAR 11 12 USE control_mod, ONLY: planet_type 12 USE strings_mod, ONLY: maxlen13 13 USE assert_eq_m, ONLY: assert_eq 14 14 USE comvert_mod, ONLY: pa,preff … … 39 39 CHARACTER(LEN=maxlen) :: msg, var, modname 40 40 INTEGER, PARAMETER :: length=100 41 INTEGER :: iq, fID, vID, idecal !, iml, jml, lml, nqt41 INTEGER :: iq, fID, vID, idecal, iqParent, iName, iZone, iPhase 42 42 REAL :: time, tab_cntrl(length) !--- RUN PARAMS TABLE 43 43 !------------------------------------------------------------------------------- … … 135 135 q(:,:,:,iq)=0. 136 136 !--- CRisi: for isotops, theoretical initialization using very simplified 137 ! Rayleigh distillation las. 138 IF(ok_isotopes.AND.iso_num(iq)>0) THEN 139 IF(zone_num(iq)==0) q(:,:,:,iq)=q(:,:,:,iqpere(iq))*tnat(iso_num(iq)) & 140 & *(q(:,:,:,iqpere(iq))/30.e-3)**(alpha_ideal(iso_num(iq))-1) 141 IF(zone_num(iq)==1) q(:,:,:,iq)=q(:,:,:,iqiso(iso_indnum(iq),phase_num(iq))) 137 ! Rayleigh distillation law. 138 iName = tracers(iq)%iso_iName 139 iZone = tracers(iq)%iso_iZone 140 iPhase= tracers(iq)%iso_iPhase 141 iqParent = tracers(iq)%iqParent 142 IF(ok_isotopes .AND. iName>0) THEN 143 IF(iZone==0) q(:,:,:,iq) = q(:,:,:,iqParent)*tnat(iName) & 144 *(q(:,:,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.) 145 IF(iZone==1) q(:,:,:,iq) = q(:,:,:,iqiso(iso_indnum(iq),iPhase)) 142 146 END IF 143 147 END DO -
LMDZ6/trunk/libf/dyn3d/dynredem.F90
r4046 r4050 7 7 USE IOIPSL 8 8 #endif 9 USE infotrac 9 USE infotrac, ONLY: nqtot, tracers, maxlen 10 10 USE netcdf, ONLY: NF90_CREATE, NF90_DEF_DIM, NF90_INQ_VARID, NF90_GLOBAL, & 11 11 NF90_CLOSE, NF90_PUT_ATT, NF90_UNLIMITED, NF90_CLOBBER, & 12 12 NF90_64BIT_OFFSET 13 13 USE dynredem_mod, ONLY: cre_var, put_var1, put_var2, err, modname, fil 14 USE strings_mod, ONLY: maxlen15 14 USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt,pa,preff, & 16 15 nivsig,nivsigs … … 167 166 ! Purpose: Write the NetCDF restart file (append). 168 167 !------------------------------------------------------------------------------- 169 USE infotrac 168 USE infotrac, ONLY: nqtot, tracers, type_trac, maxlen 170 169 USE control_mod 171 170 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_GET_VAR, NF90_INQ_VARID, & -
LMDZ6/trunk/libf/dyn3d/iniacademic.F90
r3976 r4050 5 5 6 6 USE filtreg_mod, ONLY: inifilr 7 USE infotrac 7 USE infotrac, ONLY: nqtot,niso_possibles,ok_isotopes,ok_iso_verif,tnat,alpha_ideal, & 8 & iqiso,iso_indnum, tracers 8 9 USE control_mod, ONLY: day_step,planet_type 9 10 #ifdef CPP_IOIPSL … … 61 62 real tetastrat ! potential temperature in the stratosphere, in K 62 63 real tetajl(jjp1,llm) 63 INTEGER i,j,l,lsup,ij 64 INTEGER i,j,l,lsup,ij, iq 64 65 65 66 REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T … … 70 71 71 72 real zz,ran1 72 integer idum 73 integer idum, iqParent, iName, iZone, iPhase 73 74 74 75 REAL zdtvr … … 275 276 if (planet_type=="earth") then 276 277 ! Earth: first two tracers will be water 277 do i =1,nqtot278 if (i == 1) q(:,:,i)=1.e-10279 if (i == 2) q(:,:,i)=1.e-15280 if (i.gt.2) q(:,:,i)=0.278 do iq=1,nqtot 279 q(:,:,iq)=0. 280 IF(tracers(iq)%name == 'H2Ov') q(:,:,iq)=1.e-10 281 IF(tracers(iq)%name == 'H2Ol') q(:,:,iq)=1.e-15 281 282 282 283 ! CRisi: init des isotopes 283 284 ! distill de Rayleigh très simplifiée 284 if (ok_isotopes) then 285 if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then 286 q(:,:,i)=q(:,:,iqpere(i)) & 287 & *tnat(iso_num(i)) & 288 & *(q(:,:,iqpere(i))/30.e-3) & 289 & **(alpha_ideal(iso_num(i))-1) 290 endif 291 if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then 292 q(:,:,i)=q(:,:,iqiso(iso_indnum(i),phase_num(i))) 293 endif 285 iName = tracers(iq)%iso_iName 286 iZone = tracers(iq)%iso_iZone 287 iPhase= tracers(iq)%iso_iPhase 288 iqParent = tracers(iq)%iqParent 289 if (ok_isotopes .AND. iName > 0) then 290 if (iZone == 0) q(:,:,iq) = q(:,:,iqParent)*tnat(iName) & 291 *(q(:,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1) 292 if (iZone == 1) q(:,:,iq) = q(:,:,iqiso(iso_indnum(iq),iPhase)) 294 293 endif !if (ok_isotopes) then 295 294 -
LMDZ6/trunk/libf/dyn3d/vlsplt.F
r4008 r4050 4 4 5 5 SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq) 6 USE infotrac, ONLY: nqtot, nqdesc,iqfils6 USE infotrac, ONLY: nqtot,tracers 7 7 c 8 8 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 83 83 CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1) 84 84 85 if (nqdesc(iq).gt.0) then 86 do ifils=1,nqdesc(iq) 87 iq2=iqfils(ifils,iq) 88 CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1) 89 enddo 90 endif !if (nqfils(iq).gt.0) then 85 do ifils=1,tracers(iq)%nqDescen 86 iq2=tracers(iq)%iqDescen(ifils) 87 CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1) 88 enddo 91 89 92 90 cprint*,'Entree vlx1' … … 122 120 ENDDO 123 121 ! CRisi: aussi pour les fils 124 if (nqdesc(iq).gt.0) then 125 do ifils=1,nqdesc(iq) 126 iq2=iqfils(ifils,iq) 122 do ifils=1,tracers(iq)%nqDescen 123 iq2=tracers(iq)%iqDescen(ifils) 127 124 DO l=1,llm 128 DO ij=1,ip1jmp1129 q(ij,l,iq2)=zq(ij,l,iq2)130 ENDDO131 DO ij=1,ip1jm+1,iip1125 DO ij=1,ip1jmp1 126 q(ij,l,iq2)=zq(ij,l,iq2) 127 ENDDO 128 DO ij=1,ip1jm+1,iip1 132 129 q(ij+iim,l,iq2)=q(ij,l,iq2) 133 ENDDO130 ENDDO 134 131 ENDDO 135 enddo !do ifils=1,nqdesc(iq) 136 endif ! if (nqdesc(iq).gt.0) then 132 enddo 137 133 138 134 RETURN 139 135 END 140 136 RECURSIVE SUBROUTINE vlx(q,pente_max,masse,u_m,iq) 141 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils, ! CRisi137 USE infotrac, ONLY : nqtot,tracers, ! CRisi 142 138 & qperemin,masseqmin,ratiomin ! MVals et CRisi 143 139 … … 449 445 ! CRisi: appel récursif de l'advection sur les fils. 450 446 ! Il faut faire ça avant d'avoir mis à jour q et masse 451 !write(*,*) 'vlsplt 326: iq,nq fils(iq)=',iq,nqfils(iq)447 !write(*,*) 'vlsplt 326: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen 452 448 453 if (nqdesc(iq).gt.0) then 454 do ifils=1,nqdesc(iq) 455 iq2=iqfils(ifils,iq) 456 DO l=1,llm 449 do ifils=1,tracers(iq)%nqDescen 450 iq2=tracers(iq)%iqDescen(ifils) 451 DO l=1,llm 457 452 DO ij=iip2,ip1jm 458 ! On a besoin de q et masse seulement entre iip2 et ip1jm459 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)460 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)461 !Mvals: veiller a ce qu'on n'ait pas de denominateur nul462 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)463 if (q(ij,l,iq).gt.qperemin) then464 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)465 else466 Ratio(ij,l,iq2)=ratiomin467 endif453 ! On a besoin de q et masse seulement entre iip2 et ip1jm 454 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 455 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 456 !Mvals: veiller a ce qu'on n'ait pas de denominateur nul 457 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin) 458 if (q(ij,l,iq).gt.qperemin) then 459 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 460 else 461 Ratio(ij,l,iq2)=ratiomin 462 endif 468 463 enddo 469 enddo 470 enddo !do ifils=1,nqdesc(iq) 471 do ifils=1,nqfils(iq) 472 iq2=iqfils(ifils,iq) 473 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 474 enddo !do ifils=1,nqfils(iq) 475 endif !if (nqfils(iq).gt.0) then 464 enddo 465 enddo 466 do ifils=1,tracers(iq)%nqDescen 467 iq2=tracers(iq)%iqDescen(ifils) 468 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 469 enddo 476 470 ! end CRisi 477 471 … … 498 492 ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio 499 493 ! puis on boucle en longitude 500 if (nqfils(iq).gt.0) then 501 do ifils=1,nqdesc(iq) 502 iq2=iqfils(ifils,iq) 503 DO l=1,llm 494 do ifils=1,tracers(iq)%nqDescen 495 iq2=tracers(iq)%iqDescen(ifils) 496 DO l=1,llm 504 497 DO ij=iip2+1,ip1jm 505 498 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) … … 508 501 q(ij-iim,l,iq2)=q(ij,l,iq2) 509 502 enddo ! DO ij=ijb+iip1-1,ije,iip1 510 enddo !DO l=1,llm 511 enddo !do ifils=1,nqdesc(iq) 512 endif !if (nqfils(iq).gt.0) then 503 enddo !DO l=1,llm 504 enddo 513 505 514 506 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) … … 519 511 END 520 512 RECURSIVE SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq) 521 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils, ! CRisi513 USE infotrac, ONLY : nqtot,tracers, ! CRisi 522 514 & qperemin,masseqmin,ratiomin ! MVals et CRisi 523 515 c … … 778 770 ! CRisi: appel récursif de l'advection sur les fils. 779 771 ! Il faut faire ça avant d'avoir mis à jour q et masse 780 !write(*,*) 'vly 689: iq,nq fils(iq)=',iq,nqfils(iq)772 !write(*,*) 'vly 689: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen 781 773 782 if (nqfils(iq).gt.0) then 783 do ifils=1,nqdesc(iq) 784 iq2=iqfils(ifils,iq) 785 DO l=1,llm 786 DO ij=1,ip1jmp1 787 ! attention, chaque fils doit avoir son masseq, sinon, le 1er 788 ! fils ecrase le masseq de ses freres. 789 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 790 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 791 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 792 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin) 793 if (q(ij,l,iq).gt.qperemin) then 794 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 795 else 796 Ratio(ij,l,iq2)=ratiomin 797 endif 774 do ifils=1,tracers(iq)%nqDescen 775 iq2=tracers(iq)%iqDescen(ifils) 776 DO l=1,llm 777 DO ij=1,ip1jmp1 778 ! attention, chaque fils doit avoir son masseq, sinon, le 1er 779 ! fils ecrase le masseq de ses freres. 780 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 781 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 782 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 783 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin) 784 if (q(ij,l,iq).gt.qperemin) then 785 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 786 else 787 Ratio(ij,l,iq2)=ratiomin 788 endif 798 789 enddo 799 enddo 800 enddo !do ifils=1,nqdesc(iq) 801 802 do ifils=1,nqfils(iq) 803 iq2=iqfils(ifils,iq) 804 call vly(Ratio,pente_max,masseq,qbyv,iq2) 805 enddo !do ifils=1,nqfils(iq) 806 endif !if (nqfils(iq).gt.0) then 790 enddo 791 enddo 792 793 do ifils=1,tracers(iq)%nqDescen 794 iq2=tracers(iq)%iqDescen(ifils) 795 call vly(Ratio,pente_max,masseq,qbyv,iq2) 796 enddo 807 797 808 798 DO l=1,llm … … 872 862 873 863 ! retablir les fils en rapport de melange par rapport a l'air: 874 if (nqfils(iq).gt.0) then 875 do ifils=1,nqdesc(iq) 876 iq2=iqfils(ifils,iq) 877 DO l=1,llm 864 do ifils=1,tracers(iq)%nqDescen 865 iq2=tracers(iq)%iqDescen(ifils) 866 DO l=1,llm 878 867 DO ij=1,ip1jmp1 879 868 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 880 869 enddo 881 enddo 882 enddo !do ifils=1,nqdesc(iq) 883 endif !if (nqfils(iq).gt.0) then 870 enddo 871 enddo 884 872 885 873 !write(*,*) 'vly 853: sortie' … … 888 876 END 889 877 RECURSIVE SUBROUTINE vlz(q,pente_max,masse,w,iq) 890 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils, ! CRisi878 USE infotrac, ONLY : nqtot,tracers, ! CRisi 891 879 & qperemin,masseqmin,ratiomin ! MVals et CRisi 892 880 c … … 1009 997 ! CRisi: appel récursif de l'advection sur les fils. 1010 998 ! Il faut faire ça avant d'avoir mis à jour q et masse 1011 !write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq) 1012 if (nqfils(iq).gt.0) then 1013 do ifils=1,nqdesc(iq) 1014 iq2=iqfils(ifils,iq) 1015 DO l=1,llm 999 !write(*,*) 'vlsplt 942: iq,nqChilds(iq)=',iq,nqChilds(iq) 1000 do ifils=1,tracers(iq)%nqDescen 1001 iq2=tracers(iq)%iqDescen(ifils) 1002 DO l=1,llm 1016 1003 DO ij=1,ip1jmp1 1017 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)1018 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)1019 !MVals: veiller a ce qu'on n'ait pas de denominateur nul1020 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)1021 if (q(ij,l,iq).gt.qperemin) then1022 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)1023 else1024 Ratio(ij,l,iq2)=ratiomin1025 endif1004 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 1005 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 1006 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 1007 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin) 1008 if (q(ij,l,iq).gt.qperemin) then 1009 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 1010 else 1011 Ratio(ij,l,iq2)=ratiomin 1012 endif 1026 1013 enddo 1027 1028 enddo !do ifils=1,nqdesc(iq)1014 enddo 1015 enddo 1029 1016 1030 do ifils=1,nqfils(iq) 1031 iq2=iqfils(ifils,iq) 1032 call vlz(Ratio,pente_max,masseq,wq,iq2) 1033 enddo !do ifils=1,nqfils(iq) 1034 endif !if (nqfils(iq).gt.0) then 1017 do ifils=1,tracers(iq)%nqDescen 1018 iq2=tracers(iq)%iqDescen(ifils) 1019 call vlz(Ratio,pente_max,masseq,wq,iq2) 1020 enddo 1035 1021 ! end CRisi 1036 1022 … … 1045 1031 1046 1032 ! retablir les fils en rapport de melange par rapport a l'air: 1047 if (nqfils(iq).gt.0) then 1048 do ifils=1,nqdesc(iq) 1049 iq2=iqfils(ifils,iq) 1050 DO l=1,llm 1033 do ifils=1,tracers(iq)%nqDescen 1034 iq2=tracers(iq)%iqDescen(ifils) 1035 DO l=1,llm 1051 1036 DO ij=1,ip1jmp1 1052 1037 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 1053 1038 enddo 1054 enddo 1055 enddo !do ifils=1,nqdesc(iq) 1056 endif !if (nqfils(iq).gt.0) then 1039 enddo 1040 enddo 1057 1041 !write(*,*) 'vlsplt 1032' 1058 1042 -
LMDZ6/trunk/libf/dyn3d/vlspltqs.F
r2603 r4050 4 4 SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt, 5 5 , p,pk,teta,iq ) 6 USE infotrac, ONLY: nqtot, nqdesc,iqfils6 USE infotrac, ONLY: nqtot,tracers 7 7 c 8 8 c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron … … 121 121 CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1) 122 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) 123 do ifils=1,tracers(iq)%nqDescen 124 iq2=tracers(iq)%iqDescen(ifils) 126 125 CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1) 127 enddo 128 endif !if (nqfils(iq).gt.0) then 126 enddo 129 127 130 128 c call minmaxq(zq,qmin,qmax,'avant vlxqs ') … … 162 160 ENDDO 163 161 ! CRisi: aussi pour les fils 164 if (nqdesc(iq).gt.0) then 165 do ifils=1,nqdesc(iq) 166 iq2=iqfils(ifils,iq) 162 do ifils=1,tracers(iq)%nqDescen 163 iq2=tracers(iq)%iqDescen(ifils) 167 164 DO l=1,llm 168 DO ij=1,ip1jmp1169 q(ij,l,iq2)=zq(ij,l,iq2)170 ENDDO171 DO ij=1,ip1jm+1,iip1165 DO ij=1,ip1jmp1 166 q(ij,l,iq2)=zq(ij,l,iq2) 167 ENDDO 168 DO ij=1,ip1jm+1,iip1 172 169 q(ij+iim,l,iq2)=q(ij,l,iq2) 173 ENDDO170 ENDDO 174 171 ENDDO 175 enddo !do ifils=1,nqdesc(iq) 176 endif ! if (nqfils(iq).gt.0) then 172 enddo 177 173 !write(*,*) 'vlspltqs 183: fin de la routine' 178 174 … … 180 176 END 181 177 SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq) 182 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils ! CRisi178 USE infotrac, ONLY : nqtot,tracers ! CRisi 183 179 184 180 c … … 483 479 ! CRisi: appel récursif de l'advection sur les fils. 484 480 ! Il faut faire ça avant d'avoir mis à jour q et masse 485 !write(*,*) 'vlspltqs 326: iq,nq fils(iq)=',iq,nqfils(iq)481 !write(*,*) 'vlspltqs 326: iq,nqChilds(iq)=',iq,tracers(iq)%nqChilds 486 482 487 if (nqfils(iq).gt.0) then 488 do ifils=1,nqdesc(iq) 489 iq2=iqfils(ifils,iq) 490 DO l=1,llm 483 do ifils=1,tracers(iq)%nqDescen 484 iq2=tracers(iq)%iqDescen(ifils) 485 DO l=1,llm 491 486 DO ij=iip2,ip1jm 492 493 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)494 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)487 ! On a besoin de q et masse seulement entre iip2 et ip1jm 488 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 489 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 495 490 enddo 496 enddo 497 enddo !do ifils=1,nqdesc(iq) 498 do ifils=1,nqfils(iq) 499 iq2=iqfils(ifils,iq) 500 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 501 enddo !do ifils=1,nqfils(iq) 502 endif !if (nqfils(iq).gt.0) then 491 enddo 492 enddo 493 do ifils=1,tracers(iq)%nqChilds 494 iq2=tracers(iq)%iqDescen(ifils) 495 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 496 enddo 503 497 ! end CRisi 504 498 … … 523 517 ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio 524 518 ! puis on boucle en longitude 525 if (nqdesc(iq).gt.0) then 526 do ifils=1,nqdesc(iq) 527 iq2=iqfils(ifils,iq) 528 DO l=1,llm 519 do ifils=1,tracers(iq)%nqDescen 520 iq2=tracers(iq)%iqDescen(ifils) 521 DO l=1,llm 529 522 DO ij=iip2+1,ip1jm 530 523 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 531 524 enddo 532 525 DO ij=iip1+iip1,ip1jm,iip1 533 q(ij-iim,l,iq2)=q(ij,l,iq2) 534 enddo ! DO ij=ijb+iip1-1,ije,iip1 535 enddo !DO l=1,llm 536 enddo !do ifils=1,nqdesc(iq) 537 endif !if (nqfils(iq).gt.0) then 526 q(ij-iim,l,iq2)=q(ij,l,iq2) 527 enddo 528 enddo 529 enddo 538 530 539 531 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) … … 544 536 END 545 537 SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq) 546 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils ! CRisi538 USE infotrac, ONLY : nqtot,tracers ! CRisi 547 539 c 548 540 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 794 786 ! CRisi: appel récursif de l'advection sur les fils. 795 787 ! Il faut faire ça avant d'avoir mis à jour q et masse 796 !write(*,*) 'vlyqs 689: iq,nq fils(iq)=',iq,nqfils(iq)788 !write(*,*) 'vlyqs 689: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen 797 789 798 if (nqfils(iq).gt.0) then 799 do ifils=1,nqdesc(iq) 800 iq2=iqfils(ifils,iq) 801 DO l=1,llm 802 DO ij=1,ip1jmp1 803 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 804 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 790 do ifils=1,tracers(iq)%nqDescen 791 iq2=tracers(iq)%iqDescen(ifils) 792 DO l=1,llm 793 DO ij=1,ip1jmp1 794 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 795 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 805 796 enddo 806 enddo 807 enddo !do ifils=1,nqdesc(iq) 808 809 do ifils=1,nqfils(iq) 810 iq2=iqfils(ifils,iq) 811 !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2 812 call vly(Ratio,pente_max,masseq,qbyv,iq2) 813 enddo !do ifils=1,nqfils(iq) 814 endif !if (nqfils(iq).gt.0) then 797 enddo 798 enddo 799 do ifils=1,tracers(iq)%nqChilds 800 iq2=tracers(iq)%iqDescen(ifils) 801 !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2 802 call vly(Ratio,pente_max,masseq,qbyv,iq2) 803 enddo 815 804 816 805 DO l=1,llm … … 868 857 869 858 ! retablir les fils en rapport de melange par rapport a l'air: 870 if (nqdesc(iq).gt.0) then 871 do ifils=1,nqdesc(iq) 872 iq2=iqfils(ifils,iq) 873 DO l=1,llm 859 do ifils=1,tracers(iq)%nqDescen 860 iq2=tracers(iq)%iqDescen(ifils) 861 DO l=1,llm 874 862 DO ij=1,ip1jmp1 875 863 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 876 864 enddo 877 enddo 878 enddo !do ifils=1,nqdesc(iq) 879 endif !if (nqfils(iq).gt.0) then 865 enddo 866 enddo 880 867 !write(*,*) 'vly 879' 881 868 -
LMDZ6/trunk/libf/dyn3d_common/infotrac.F90
r4046 r4050 32 32 TYPE(isot_type), TARGET, SAVE, ALLOCATABLE :: isotopes(:) !=== ISOTOPES PARAMETERS VECTOR 33 33 34 ! iadv : index of trasport schema for each tracer35 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iadv36 37 34 ! niadv : vector keeping the coorspondance between all tracers(nqtot) treated in the 38 35 ! dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code. 39 36 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: niadv ! equivalent dyn / physique 40 37 41 ! CRisi: tableaux de fils42 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqfils43 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqdesc ! nombres de fils + nombre de tous les petits fils sur toutes les generations44 INTEGER, SAVE :: nqdesc_tot45 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqfils46 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iqpere47 38 REAL :: qperemin,masseqmin,ratiomin ! MVals et CRisi 48 39 PARAMETER (qperemin=1e-30,masseqmin=1e-18,ratiomin=1e-16) ! MVals … … 63 54 LOGICAL, DIMENSION(niso_possibles),SAVE :: use_iso 64 55 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqiso ! donne indice iq en fn de (ixt,phase) 65 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_num ! donne numero iso entre 1 et niso_possibles en fn de nqtot66 56 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_indnum ! donne numero iso entre 1 et niso effectif en fn de nqtot 67 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: zone_num ! donne numero de la zone de tracage en fn de nqtot68 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: phase_num ! donne numero de la zone de tracage en fn de nqtot69 57 INTEGER, DIMENSION(niso_possibles), SAVE :: indnum_fn_num ! donne indice entre entre 1 et niso en fonction du numero d isotope entre 1 et niso_possibles 70 58 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: index_trac ! numero ixt en fn izone, indnum entre 1 et niso … … 128 116 129 117 INTEGER :: nqtrue ! number of tracers read from tracer.def, without higer order of moment 130 INTEGER :: iq, new_iq, iiq, jq, ierr,itr 131 INTEGER :: ifils,ipere ,generation! CRisi118 INTEGER :: iq, new_iq, iiq, jq, ierr,itr, iadv 119 INTEGER :: ifils,ipere ! CRisi 132 120 LOGICAL :: continu,nouveau_traceurdef 133 121 INTEGER :: IOstatus ! gestion de la retrocompatibilite de traceur.def 134 122 CHARACTER(len=maxlen) :: tchaine 123 INTEGER, ALLOCATABLE :: iqfils(:,:) 135 124 136 125 character(len=*),parameter :: modname="infotrac_init" … … 565 554 ! 566 555 ALLOCATE(tracers(nqtot)) 567 ALLOCATE( iadv(nqtot),niadv(nqtot))556 ALLOCATE(niadv(nqtot)) 568 557 569 558 !----------------------------------------------------------------------- … … 578 567 ! Verify choice of advection schema 579 568 IF (hadv(iq)==vadv(iq)) THEN 580 iadv(new_iq)=hadv(iq)569 tracers(new_iq)%iadv=hadv(iq) 581 570 ELSE IF (hadv(iq)==10 .AND. vadv(iq)==16) THEN 582 iadv(new_iq)=11571 tracers(new_iq)%iadv=11 583 572 ELSE 584 573 WRITE(lunout,*)trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq) … … 588 577 589 578 str1=tnom_0(iq) 590 tracers(new_iq)%name=TRIM(tnom_0(iq)) 591 IF (iadv(new_iq)==0) THEN 579 tracers(new_iq)%name = TRIM(tnom_0(iq)) 580 tracers(new_iq)%parent = TRIM(tnom_transp(iq)) 581 IF (tracers(new_iq)%iadv==0) THEN 592 582 tracers(new_iq)%longName=trim(str1) 593 583 ELSE 594 tracers(new_iq)%longName=trim(tnom_0(iq))//descrq( iadv(new_iq))584 tracers(new_iq)%longName=trim(tnom_0(iq))//descrq(tracers(new_iq)%iadv) 595 585 ENDIF 596 586 597 587 ! schemas tenant compte des moments d'ordre superieur 598 588 str2=TRIM(tracers(new_iq)%longName) 599 IF ( iadv(new_iq)==20) THEN589 IF (tracers(new_iq)%iadv==20) THEN 600 590 DO jq=1,3 601 591 new_iq=new_iq+1 602 iadv(new_iq)=-20592 tracers(new_iq)%iadv=-20 603 593 tracers(new_iq)%longName=trim(str2)//txts(jq) 604 594 tracers(new_iq)%name=trim(str1)//txts(jq) 605 595 END DO 606 ELSE IF ( iadv(new_iq)==30) THEN596 ELSE IF (tracers(new_iq)%iadv==30) THEN 607 597 DO jq=1,9 608 598 new_iq=new_iq+1 609 iadv(new_iq)=-30599 tracers(new_iq)%iadv=-30 610 600 tracers(new_iq)%longName=trim(str2)//txtp(jq) 611 601 tracers(new_iq)%name=trim(str1)//txtp(jq) … … 620 610 iiq=0 621 611 DO iq=1,nqtot 622 IF( iadv(iq).GE.0) THEN612 IF(tracers(iq)%iadv.GE.0) THEN 623 613 ! True tracer 624 614 iiq=iiq+1 … … 632 622 633 623 DO iq=1,nqtot 634 WRITE(lunout,*) iadv(iq),niadv(iq), ' ',trim(tracers(iq)%name),' ',trim(tracers(iq)%longName)624 WRITE(lunout,*) tracers(iq)%iadv,niadv(iq), ' ',trim(tracers(iq)%name),' ',trim(tracers(iq)%longName) 635 625 END DO 636 626 … … 640 630 ! 641 631 DO iq=1,nqtot 642 IF (iadv(iq)/=10 .AND. iadv(iq)/=14 .AND. iadv(iq)/=0) THEN 643 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ' 632 iadv=tracers(iq)%iadv 633 IF (ALL([10, 14, 0]/=iadv)) THEN 634 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv,' is not tested in this version of LMDZ' 644 635 CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1) 645 ELSE IF (iadv (iq)==14 .AND. iq/=1) THEN646 WRITE(lunout,*)trim(modname),'STOP : The option iadv=',iadv (iq),' is not tested in this version of LMDZ'636 ELSE IF (iadv==14 .AND. iq/=1) THEN 637 WRITE(lunout,*)trim(modname),'STOP : The option iadv=',iadv,' is not tested in this version of LMDZ' 647 638 CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1) 648 639 ENDIF … … 653 644 ! initialiser tous les tableaux d'indices lies aux traceurs familiaux 654 645 ! + verifier que tous les peres sont ecrits en premieres positions 655 ALLOCATE(nqfils(nqtot),nqdesc(nqtot))656 646 ALLOCATE(iqfils(nqtot,nqtot)) 657 ALLOCATE(iqpere(nqtot))658 647 nqperes=0 659 nqfils(:)=0660 nqdesc(:)=0661 648 iqfils(:,:)=0 662 iqpere(:)=0 663 nqdesc_tot=0 649 tracers(:)%iqParent=0 664 650 DO iq=1,nqtot 665 651 if (tnom_transp(iq) == 'air') then … … 667 653 WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un pere' 668 654 nqperes=nqperes+1 669 iqpere(iq)=0655 tracers(iq)%iqParent=0 670 656 else !if (tnom_transp(iq) == 'air') then 671 657 ! ceci est un fils. Qui est son pere? … … 681 667 CALL abort_gcm('infotrac_init','Un fils est son propre pere',1) 682 668 endif 683 nqfils(ipere)=nqfils(ipere)+1684 iqfils( nqfils(ipere),ipere)=iq685 iqpere(iq)=ipere669 tracers(ipere)%nqChilds = tracers(ipere)%nqChilds+1 670 iqfils(tracers(ipere)%nqChilds,ipere)=iq 671 tracers(iq)%iqParent=ipere 686 672 continu=.false. 687 673 else !if (tnom_transp(iq) == tnom_0(ipere)) then … … 697 683 enddo !DO iq=1,nqtot 698 684 WRITE(lunout,*) 'infotrac: nqperes=',nqperes 699 WRITE(lunout,*) 'nq fils=',nqfils700 WRITE(lunout,*) 'iq pere=',iqpere685 WRITE(lunout,*) 'nqChilds=',tracers(:)%nqChilds 686 WRITE(lunout,*) 'iqParent=',tracers(:)%iqParent 701 687 WRITE(lunout,*) 'iqfils=',iqfils 702 688 703 689 ! Calculer le nombre de descendants a partir de iqfils et de nbfils 704 690 DO iq=1,nqtot 705 generation=0691 tracers(iq)%iGeneration=0 706 692 continu=.true. 707 693 ifils=iq 708 694 do while (continu) 709 ipere= iqpere(ifils)695 ipere=tracers(ifils)%iqParent 710 696 if (ipere.gt.0) then 711 nqdesc(ipere)=nqdesc(ipere)+1 712 nqdesc_tot=nqdesc_tot+1 713 iqfils(nqdesc(ipere),ipere)=iq 697 tracers(ipere)%nqDescen = tracers(ipere)%nqDescen+1 698 iqfils(tracers(ipere)%nqDescen,ipere)=iq 714 699 ifils=ipere 715 generation=generation+1700 tracers(iq)%iGeneration=tracers(iq)%iGeneration+1 716 701 else !if (ipere.gt.0) then 717 702 continu=.false. 718 703 endif !if (ipere.gt.0) then 719 704 enddo !do while (continu) 720 WRITE(lunout,*) 'Le traceur ',iq,', appele ',trim(tnom_0(iq)),' est un traceur de generation: ',generation 705 WRITE(lunout,*) 'Le traceur ',iq,', appele ',trim(tnom_0(iq)), & 706 ' est un traceur de generation: ',tracers(iq)%iGeneration 721 707 enddo !DO iq=1,nqtot 722 WRITE(lunout,*) 'infotrac: nqdesc=',nqdesc 708 DO iq=1,nqtot 709 ALLOCATE(tracers(iq)%iqDescen(tracers(iq)%nqDescen)) 710 tracers(iq)%iqDescen(:) = iqfils(1:tracers(iq)%nqDescen,iq) 711 END DO 712 713 WRITE(lunout,*) 'infotrac: nqDescen=',tracers(iq)%nqDescen 723 714 WRITE(lunout,*) 'iqfils=',iqfils 724 WRITE(lunout,*) 'nq desc_tot=',nqdesc_tot715 WRITE(lunout,*) 'nqDescen_tot=',SUM(tracers(:)%nqDescen) 725 716 726 717 ! Interdire autres schemas que 10 pour les traceurs fils, et autres schemas 727 718 ! que 10 et 14 si des peres ont des fils 728 719 do iq=1,nqtot 729 if ( iqpere(iq).gt.0) then720 if (tracers(iq)%iqParent > 0) then 730 721 ! ce traceur a un pere qui n'est pas l'air 731 722 ! Seul le schema 10 est autorise 732 if (iadv(iq)/=10) then 733 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for sons' 723 iadv=tracers(iq)%iadv 724 if (iadv/=10) then 725 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv,' is not implemented for sons' 734 726 CALL abort_gcm('infotrac_init','Sons should be advected by scheme 10',1) 735 727 endif 736 728 ! Le traceur pere ne peut etre advecte que par schema 10 ou 14: 737 IF ( iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN738 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv (iq),' is not implemented for fathers'729 IF (ALL([10,14]/=tracers(tracers(iq)%iqParent)%iadv)) THEN 730 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv,' is not implemented for fathers' 739 731 CALL abort_gcm('infotrac_init','Fathers should be advected by scheme 10 ou 14',1) 740 endif !IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN741 endif !if (iqpere(iq).gt.0) the732 endif 733 endif 742 734 enddo !do iq=1,nqtot 743 735 tracers(:)%gen0Name = ancestor(tracers) !--- Name of the first generation ancestor 744 736 745 737 … … 763 755 itr=0 764 756 do iq=nqo+1, nqtot 765 if ( iso_num(iq).eq.0) then757 if (tracers(iq)%iso_iName.eq.0) then 766 758 itr=itr+1 767 759 write(*,*) 'itr=',itr 768 760 itr_indice(itr)=iq 769 endif !if ( iso_num(iq).eq.0) then761 endif !if (tracers(iq)%iso_iName.eq.0) then 770 762 enddo 771 763 if (itr.ne.nqtottr) then … … 811 803 ALLOCATE(nb_isoind(nqo)) 812 804 ALLOCATE(nb_traciso(niso_possibles,nqo)) 813 ALLOCATE(iso_num(nqtot))814 805 ALLOCATE(iso_indnum(nqtot)) 815 ALLOCATE(zone_num(nqtot))816 ALLOCATE(phase_num(nqtot))817 806 818 iso_num(:)=0819 807 iso_indnum(:)=0 820 zone_num(:)=0821 phase_num(:)=0822 808 indnum_fn_num(:)=0 823 809 use_iso(:)=.false. … … 841 827 nb_iso(ixt,phase)=nb_iso(ixt,phase)+1 842 828 nb_isoind(phase)=nb_isoind(phase)+1 843 iso_num(iq)=ixt829 tracers(iq)%iso_iName=ixt 844 830 iso_indnum(iq)=nb_isoind(phase) 845 831 indnum_fn_num(ixt)=iso_indnum(iq) 846 phase_num(iq)=phase 847 ! write(lunout,*) 'iso_num(iq)=',iso_num(iq) 848 ! write(lunout,*) 'iso_indnum(iq)=',iso_indnum(iq) 849 ! write(lunout,*) 'indnum_fn_num(ixt)=',indnum_fn_num(ixt) 850 ! write(lunout,*) 'phase_num(iq)=',phase_num(iq) 832 tracers(iq)%iso_iPhase=phase 851 833 goto 20 852 else if ( iqpere(iq).gt.0) then853 if (tnom_0( iqpere(iq)) == tnom_trac) then834 else if ( tracers(iq)%iqParent> 0) then 835 if (tnom_0(tracers(iq)%iqParent) == tnom_trac) then 854 836 ! write(lunout,*) 'Ce traceur est le fils d''un isotope' 855 837 ! c'est un traceur d'isotope 856 838 nb_traciso(ixt,phase)=nb_traciso(ixt,phase)+1 857 iso_num(iq)=ixt839 tracers(iq)%iso_iName=ixt 858 840 iso_indnum(iq)=indnum_fn_num(ixt) 859 zone_num(iq)=nb_traciso(ixt,phase) 860 phase_num(iq)=phase 861 ! write(lunout,*) 'iso_num(iq)=',iso_num(iq) 862 ! write(lunout,*) 'phase_num(iq)=',phase_num(iq) 863 ! write(lunout,*) 'zone_num(iq)=',zone_num(iq) 841 tracers(iq)%iso_iZone=nb_traciso(ixt,phase) 842 tracers(iq)%iso_iPhase=phase 864 843 goto 20 865 endif !if (tnom_0( iqpere(iq)) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then844 endif !if (tnom_0(tracers(iq)%iqParent) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then 866 845 endif !IF (tnom_0(iq) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then 867 846 enddo !do ixt= niso_possibles … … 869 848 20 continue 870 849 enddo !do iq=1,nqtot 871 872 ! write(lunout,*) 'iso_num=',iso_num873 ! write(lunout,*) 'iso_indnum=',iso_indnum874 ! write(lunout,*) 'zone_num=',zone_num875 ! write(lunout,*) 'phase_num=',phase_num876 ! write(lunout,*) 'indnum_fn_num=',indnum_fn_num877 850 878 851 do ixt= 1,niso_possibles … … 926 899 927 900 ! flags isotopiques: 928 if (niso.gt.0) then 929 ok_isotopes=.true. 930 else 931 ok_isotopes=.false. 932 endif 901 ok_isotopes = niso > 0 933 902 ! WRITE(lunout,*) 'ok_isotopes=',ok_isotopes 934 903 … … 955 924 iqiso(:,:)=0 956 925 do iq=1,nqtot 957 if ( iso_num(iq).gt.0) then958 ixt=iso_indnum(iq)+ zone_num(iq)*niso959 iqiso(ixt, phase_num(iq))=iq926 if (tracers(iq)%iso_iName > 0) then 927 ixt=iso_indnum(iq)+tracers(iq)%iso_iZone*niso 928 iqiso(ixt,tracers(iq)%iso_iPhase)=iq 960 929 endif 961 930 enddo -
LMDZ6/trunk/libf/dyn3d_common/iso_verif_dyn.F
r2270 r4050 64 64 function iso_verif_aberrant_nostop 65 65 : (x,iso,q,err_msg) 66 USE infotrac 66 USE infotrac, ONLY: tnat 67 67 implicit none 68 68 -
LMDZ6/trunk/libf/dyn3dmem/advtrac_loc.F
r4038 r4050 24 24 USE Vampir 25 25 USE times 26 USE infotrac, ONLY: nqtot, iadv, ok_iso_verif26 USE infotrac, ONLY: nqtot, tracers, ok_iso_verif 27 27 USE control_mod, ONLY: iapp_tracvl, day_step, planet_type 28 28 USE advtrac_mod, ONLY: finmasse … … 74 74 DATA fill/.true./ 75 75 DATA dum/.true./ 76 integer ijb,ije,ijbu,ijbv,ijeu,ijev,j 76 integer ijb,ije,ijbu,ijbv,ijeu,ijev,j, iadv 77 77 type(Request),SAVE :: testRequest 78 78 !$OMP THREADPRIVATE(testRequest) … … 152 152 153 153 !write(*,*) 'advtrac 157: appel de vlspltgen_loc' 154 call vlspltgen_loc( q,iadv, 2., massem, wg , 155 * pbarug,pbarvg,dtvr,p, 154 call vlspltgen_loc( q, 2., massem, wg,pbarug,pbarvg,dtvr,p, 156 155 * pk,teta ) 157 156 … … 169 168 do iq=1,nqtot 170 169 c call clock(t_initial) 171 if(iadv(iq) == 0) cycle 170 iadv = tracers(iq)%iadv 171 SELECT CASE(iadv) 172 CASE(0); CYCLE 173 CASE(10) 172 174 c ---------------------------------------------------------------- 173 175 c Schema de Van Leer I MUSCL 174 176 c ---------------------------------------------------------------- 175 if(iadv(iq).eq.10) THEN176 177 177 178 !LF call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr) 178 179 180 CASE(14) 179 181 c ---------------------------------------------------------------- 180 182 c Schema "pseudo amont" + test sur humidite specifique 181 183 C pour la vapeur d'eau. F. Codron 182 184 c ---------------------------------------------------------------- 183 else if(iadv(iq).eq.14) then184 185 c 185 186 cym stop 'advtrac : appel à vlspltqs :schema non parallelise' 186 187 !LF CALL vlspltqs_p( q(1,1,1), 2., massem, wg , 187 188 !LF * pbarug,pbarvg,dtvr,p,pk,teta ) 189 CASE(12) 188 190 c ---------------------------------------------------------------- 189 191 c Schema de Frederic Hourdin 190 192 c ---------------------------------------------------------------- 191 else if(iadv(iq).eq.12) then192 193 stop 'advtrac : schema non parallelise' 193 194 c Pas de temps adaptatif 194 call adaptdt(iadv (iq),dtbon,n,pbarug,massem)195 call adaptdt(iadv,dtbon,n,pbarug,massem) 195 196 if (n.GT.1) then 196 197 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', … … 200 201 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1) 201 202 end do 202 else if(iadv(iq).eq.13) then203 CASE(13) 203 204 stop 'advtrac : schema non parallelise' 204 205 c Pas de temps adaptatif 205 call adaptdt(iadv (iq),dtbon,n,pbarug,massem)206 call adaptdt(iadv,dtbon,n,pbarug,massem) 206 207 if (n.GT.1) then 207 208 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', … … 211 212 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2) 212 213 end do 214 CASE(20) 213 215 c ---------------------------------------------------------------- 214 216 c Schema de pente SLOPES 215 217 c ---------------------------------------------------------------- 216 else if (iadv(iq).eq.20) then217 218 stop 'advtrac : schema non parallelise' 218 219 219 220 call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0) 220 221 222 CASE(30) 221 223 c ---------------------------------------------------------------- 222 224 c Schema de Prather 223 225 c ---------------------------------------------------------------- 224 else if (iadv(iq).eq.30) then225 226 stop 'advtrac : schema non parallelise' 226 227 c Pas de temps adaptatif 227 call adaptdt(iadv (iq),dtbon,n,pbarug,massem)228 call adaptdt(iadv,dtbon,n,pbarug,massem) 228 229 if (n.GT.1) then 229 230 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', … … 232 233 call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, 233 234 s n,dtbon) 235 CASE(11,16,17,18) 234 236 c ---------------------------------------------------------------- 235 237 c Schemas PPM Lin et Rood 236 238 c ---------------------------------------------------------------- 237 else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.238 s iadv(iq).LE.18)) then239 239 240 240 stop 'advtrac : schema non parallelise' … … 242 242 c Test sur le flux horizontal 243 243 c Pas de temps adaptatif 244 call adaptdt(iadv (iq),dtbon,n,pbarug,massem)244 call adaptdt(iadv,dtbon,n,pbarug,massem) 245 245 if (n.GT.1) then 246 246 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', … … 273 273 c VL (version PPM) horiz. et PPM vert. 274 274 c--------------------------------------------------------------------- 275 if (iadv(iq).eq.11) then 275 SELECT CASE(iadv) 276 CASE(11) 276 277 c Ss-prg PPM3d de Lin 277 278 call ppm3d(1,qppm(1,1,iq), … … 281 282 s fill,dum,220.) 282 283 284 CASE(16) 283 285 c---------------------------------------------------------------------- 284 286 c Monotonic PPM 285 287 c---------------------------------------------------------------------- 286 else if (iadv(iq).eq.16) then287 288 c Ss-prg PPM3d de Lin 288 289 call ppm3d(1,qppm(1,1,iq), … … 293 294 c--------------------------------------------------------------------- 294 295 296 CASE(17) 295 297 c--------------------------------------------------------------------- 296 298 c Semi Monotonic PPM 297 299 c--------------------------------------------------------------------- 298 else if (iadv(iq).eq.17) then299 300 c Ss-prg PPM3d de Lin 300 301 call ppm3d(1,qppm(1,1,iq), … … 305 306 c--------------------------------------------------------------------- 306 307 308 CASE(18) 307 309 c--------------------------------------------------------------------- 308 310 c Positive Definite PPM 309 311 c--------------------------------------------------------------------- 310 else if (iadv(iq).eq.18) then311 312 c Ss-prg PPM3d de Lin 312 313 call ppm3d(1,qppm(1,1,iq), … … 316 317 s fill,dum,220.) 317 318 c--------------------------------------------------------------------- 318 endif319 END SELECT 319 320 enddo 320 321 c----------------------------------------------------------------- … … 322 323 c----------------------------------------------------------------- 323 324 call interpost(q(1,1,iq),qppm(1,1,iq)) 324 endif325 END SELECT 325 326 c---------------------------------------------------------------------- 326 327 -
LMDZ6/trunk/libf/dyn3dmem/advtrac_mod.F90
r1907 r4050 9 9 USE allocate_field_mod 10 10 USE parallel_lmdz 11 USE infotrac12 11 USE vlspltgen_mod 13 12 IMPLICIT NONE -
LMDZ6/trunk/libf/dyn3dmem/caladvtrac_mod.F90
r3435 r4050 23 23 USE allocate_field_mod 24 24 USE parallel_lmdz 25 USE infotrac 25 USE infotrac, ONLY: nqtot 26 26 USE advtrac_mod, ONLY : advtrac_allocate 27 27 USE groupe_mod -
LMDZ6/trunk/libf/dyn3dmem/call_calfis_mod.F90
r3435 r4050 37 37 USE parallel_lmdz 38 38 USE dimensions_mod 39 USE infotrac 39 USE infotrac, ONLY: nqtot 40 40 IMPLICIT NONE 41 41 TYPE(distrib),POINTER :: d … … 80 80 USE Bands 81 81 USE vampir 82 USE infotrac 82 USE infotrac, ONLY: nqtot 83 83 USE control_mod 84 84 USE write_field_loc -
LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F
r4038 r4050 1 1 subroutine check_isotopes(q,ijb,ije,err_msg) 2 USE infotrac 2 USE infotrac, ONLY: nqtot, nqo, niso, ntraciso, ntraceurs_zone, 3 & ok_isotopes, ok_isotrac, use_iso, 4 & iqiso, indnum_fn_num, index_trac, tnat 3 5 USE parallel_lmdz 4 6 implicit none -
LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90
r4046 r4050 7 7 !------------------------------------------------------------------------------- 8 8 USE parallel_lmdz 9 USE infotrac 9 USE infotrac, ONLY: nqtot, tracers, iqiso, iso_indnum, tnat, alpha_ideal, ok_isotopes, maxlen 10 10 USE netcdf, ONLY: NF90_OPEN, NF90_INQUIRE_DIMENSION, NF90_INQ_VARID, & 11 11 NF90_NOWRITE, NF90_CLOSE, NF90_INQUIRE_VARIABLE, NF90_GET_VAR, NF90_NoErr 12 12 USE control_mod, ONLY: planet_type 13 USE strings_mod, ONLY: maxlen14 13 USE assert_eq_m, ONLY: assert_eq 15 14 USE comvert_mod, ONLY: pa,preff … … 42 41 CHARACTER(LEN=maxlen) :: msg, var, modname 43 42 INTEGER, PARAMETER :: length=100 44 INTEGER :: iq, fID, vID, idecal, ierr 43 INTEGER :: iq, fID, vID, idecal, ierr, iqParent, iName, iZone, iPhase 45 44 REAL :: time, tab_cntrl(length) !--- RUN PARAMS TABLE 46 45 REAL, ALLOCATABLE :: vcov_glo(:,:),masse_glo(:,:), ps_glo(:) … … 174 173 !--- CRisi: for isotops, theoretical initialization using very simplified 175 174 ! Rayleigh distillation las. 176 IF(ok_isotopes.AND.iso_num(iq)>0) THEN 177 IF(zone_num(iq)==0) q(:,:,iq)=q(:,:,iqpere(iq))*tnat(iso_num(iq)) & 178 & *(q(:,:,iqpere(iq))/30.e-3)**(alpha_ideal(iso_num(iq))-1) 179 IF(zone_num(iq)==1) q(:,:,iq)=q(:,:,iqiso(iso_indnum(iq),phase_num(iq))) 175 iName = tracers(iq)%iso_iName 176 iZone = tracers(iq)%iso_iZone 177 iPhase= tracers(iq)%iso_iPhase 178 iqParent = tracers(iq)%iqParent 179 IF(ok_isotopes .AND. iName>0) THEN 180 IF(iZone==0) q(:,:,iq) = q(:,:,iqParent)*tnat(iName) & 181 & *(q(:,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.) 182 IF(iZone==1) q(:,:,iq) = q(:,:,iqiso(iso_indnum(iq),iPhase)) 180 183 END IF 181 184 END DO -
LMDZ6/trunk/libf/dyn3dmem/dynredem_loc.F90
r4046 r4050 9 9 USE parallel_lmdz 10 10 USE mod_hallo 11 USE infotrac 11 USE infotrac, ONLY: nqtot, tracers, maxlen 12 12 USE netcdf, ONLY: NF90_CREATE, NF90_DEF_DIM, NF90_INQ_VARID, NF90_GLOBAL, & 13 13 NF90_CLOSE, NF90_PUT_ATT, NF90_UNLIMITED, NF90_CLOBBER, & 14 14 NF90_64BIT_OFFSET 15 15 USE dynredem_mod, ONLY: cre_var, put_var, err, modname, fil 16 USE strings_mod, ONLY: maxlen17 16 USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt,pa,preff, & 18 17 nivsig,nivsigs … … 175 174 USE parallel_lmdz 176 175 USE mod_hallo 177 USE infotrac 176 USE infotrac, ONLY: nqtot, tracers, type_trac, maxlen 178 177 USE control_mod 179 178 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_GET_VAR, NF90_INQ_VARID, & -
LMDZ6/trunk/libf/dyn3dmem/gcm.F90
r3579 r4050 9 9 USE mod_const_mpi, ONLY: init_const_mpi 10 10 USE parallel_lmdz 11 USE infotrac 11 USE infotrac, ONLY: nqtot, infotrac_init 12 12 !#ifdef CPP_PHYS 13 13 ! USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys -
LMDZ6/trunk/libf/dyn3dmem/groupe_mod.F90
r1907 r4050 10 10 USE allocate_field_mod 11 11 USE parallel_lmdz 12 USE infotrac12 ! USE infotrac 13 13 USE advtrac_mod, ONLY : advtrac_allocate 14 14 IMPLICIT NONE -
LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90
r3976 r4050 7 7 use exner_hyb_m, only: exner_hyb 8 8 use exner_milieu_m, only: exner_milieu 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_num9 USE infotrac, ONLY: nqtot, niso_possibles, ok_isotopes, ok_iso_verif, tnat, alpha_ideal, & 10 iqiso, tracers, iso_indnum 11 11 USE control_mod, ONLY: day_step,planet_type 12 12 USE parallel_lmdz, ONLY: ijb_u, ije_u, ijb_v, ije_v … … 67 67 real tetastrat ! potential temperature in the stratosphere, in K 68 68 real tetajl(jjp1,llm) 69 INTEGER i,j,l,lsup,ij 69 INTEGER i,j,l,lsup,ij, iq, iName, iZone, iPhase, iqParent 70 70 71 71 REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T … … 282 282 ! Earth: first two tracers will be water 283 283 284 do i =1,nqtot285 if (i == 1) q(ijb_u:ije_u,:,i)=1.e-10286 if ( i == 2) q(ijb_u:ije_u,:,i)=1.e-15287 if ( i.gt.2) q(ijb_u:ije_u,:,i)=0.284 do iq=1,nqtot 285 q(ijb_u:ije_u,:,iq)=0. 286 if (tracers(iq)%name == 'H2Ov') q(ijb_u:ije_u,:,iq)=1.e-10 287 if (tracers(iq)%name == 'H2Ol') q(ijb_u:ije_u,:,iq)=1.e-15 288 288 289 289 ! CRisi: init des isotopes 290 290 ! distill de Rayleigh très simplifiée 291 if (ok_isotopes) then 292 if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then 293 q(ijb_u:ije_u,:,i)=q(ijb_u:ije_u,:,iqpere(i)) & 294 & *tnat(iso_num(i)) & 295 & *(q(ijb_u:ije_u,:,iqpere(i))/30.e-3) & 296 & **(alpha_ideal(iso_num(i))-1) 297 endif 298 if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then 299 q(ijb_u:ije_u,:,i)=q(ijb_u:ije_u,:,iqiso(iso_indnum(i),phase_num(i))) 300 endif 291 iName = tracers(iq)%iso_iName 292 iZone = tracers(iq)%iso_iZone 293 iPhase= tracers(iq)%iso_iPhase 294 iqParent = tracers(iq)%iqParent 295 if (ok_isotopes .AND. iName > 0) then 296 if (iZone == 0) q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat(iName) & 297 *(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1) 298 if (iZone == 1) q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqiso(iso_indnum(iq),iPhase)) 301 299 endif !if (ok_isotopes) then 302 300 -
LMDZ6/trunk/libf/dyn3dmem/initdynav_loc.F
r4046 r4050 11 11 use Write_field 12 12 use misc_mod 13 USE infotrac13 ! USE infotrac 14 14 use com_io_dyn_mod, only : histaveid,histvaveid,histuaveid, & 15 15 & dynhistave_file,dynhistvave_file,dynhistuave_file -
LMDZ6/trunk/libf/dyn3dmem/inithist_loc.F
r4046 r4050 11 11 use Write_field 12 12 use misc_mod 13 USE infotrac14 13 use com_io_dyn_mod, only : histid,histvid,histuid, & 15 14 & dynhist_file,dynhistv_file,dynhistu_file -
LMDZ6/trunk/libf/dyn3dmem/vlsplt_loc.F
r3800 r4050 14 14 c -------------------------------------------------------------------- 15 15 USE parallel_lmdz 16 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils, ! CRisi &16 USE infotrac, ONLY : nqtot,tracers, ! CRisi & 17 17 & qperemin,masseqmin,ratiomin ! MVals et CRisi 18 18 IMPLICIT NONE … … 330 330 ! Il faut faire ça avant d'avoir mis à jour q et masse 331 331 332 if (nqfils(iq).gt.0) then 333 do ifils=1,nqdesc(iq) 334 !do ifils=1,nqfils(iq) ! modif C Risi 22nov2020 332 do ifils=1,tracers(iq)%nqDescen 335 333 ! attention: comme Ratio est utilisé comme q dans l'appel 336 334 ! recursif, il doit contenir à lui seul tous les indices de tous 337 335 ! les descendants! 338 iq2=iqfils(ifils,iq)339 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 340 336 iq2=tracers(iq)%iqDescen(ifils) 337 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 338 DO l=1,llm 341 339 DO ij=ijb,ije 342 ! On a besoin de q et masse seulement entre ijb et ije. On ne343 ! les calcule donc que de ijb à ije344 !MVals: veiller a ce qu'on n'ait pas de denominateur nul345 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)346 if (q(ij,l,iq).gt.qperemin) then ! modif 13 nov 2020347 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)348 else349 Ratio(ij,l,iq2)=ratiomin350 endif340 ! On a besoin de q et masse seulement entre ijb et ije. On ne 341 ! les calcule donc que de ijb à ije 342 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 343 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin) 344 if (q(ij,l,iq).gt.qperemin) then ! modif 13 nov 2020 345 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 346 else 347 Ratio(ij,l,iq2)=ratiomin 348 endif 351 349 enddo 352 enddo 353 c$OMP END DO NOWAIT 354 enddo !do ifils=1,nqdesc(iq) 355 do ifils=1,nqfils(iq) 356 iq2=iqfils(ifils,iq) 357 call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2) 358 enddo !do ifils=1,nqfils(iq) 359 endif !if (nqfils(iq).gt.0) then 350 enddo 351 c$OMP END DO NOWAIT 352 enddo !do ifils=1,tracers(iq)%nqDescen 353 do ifils=1,tracers(iq)%nqChilds 354 iq2=tracers(iq)%iqDescen(ifils) 355 call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2) 356 enddo 360 357 ! end CRisi 361 358 … … 383 380 ! On calcule q entre ijb+1 et ije -> on fait pareil pour ratio 384 381 ! puis on boucle en longitude 385 if (nqfils(iq).gt.0) then 386 do ifils=1,nqdesc(iq) 387 !do ifils=1,nqfils(iq) ! modif C Risi 22nov2020 388 iq2=iqfils(ifils,iq) 382 do ifils=1,tracers(iq)%nqDescen 383 iq2=tracers(iq)%iqDescen(ifils) 389 384 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 390 385 DO l=1,llm 391 386 DO ij=ijb+1,ije 392 387 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 393 388 enddo 394 389 DO ij=ijb+iip1-1,ije,iip1 395 q(ij-iim,l,iq2)=q(ij,l,iq2) 396 enddo ! DO ij=ijb+iip1-1,ije,iip1 397 enddo !DO l=1,llm 398 c$OMP END DO NOWAIT 399 enddo !do ifils=1,nqdesc(iq) 400 endif !if (nqfils(iq).gt.0) then 390 q(ij-iim,l,iq2)=q(ij,l,iq2) 391 enddo 392 enddo 393 c$OMP END DO NOWAIT 394 enddo 401 395 402 396 !write(*,*) 'vlsplt 399: iq,ijb_x=',iq,ijb_x … … 422 416 c -------------------------------------------------------------------- 423 417 USE parallel_lmdz 424 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils, ! CRisi &418 USE infotrac, ONLY : nqtot,tracers, ! CRisi & 425 419 & qperemin,masseqmin,ratiomin ! MVals et CRisi 426 420 USE comconst_mod, ONLY: pi … … 732 726 ! CRisi: appel récursif de l'advection sur les fils. 733 727 ! Il faut faire ça avant d'avoir mis à jour q et masse 734 !write(*,*) 'vly 689: iq,nq fils(iq)=',iq,nqfils(iq)728 !write(*,*) 'vly 689: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen 735 729 736 730 ijb=ij_begin-2*iip1 … … 743 737 if (pole_sud) ijem=ij_end 744 738 745 if (nqfils(iq).gt.0) then 746 do ifils=1,nqdesc(iq) 747 !do ifils=1,nqfils(iq) ! modif C Risi 22nov2020 748 iq2=iqfils(ifils,iq) 739 do ifils=1,tracers(iq)%nqDescen 740 iq2=tracers(iq)%iqDescen(ifils) 749 741 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 750 751 752 742 DO l=1,llm 743 ! modif des bornes: CRisi 16 nov 2020 744 ! d'abord masse avec bornes corrigées 753 745 DO ij=ijbm,ijem 754 755 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)756 enddo !DO ij=ijbm,ijem746 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 747 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin) 748 enddo 757 749 758 750 ! ensuite Ratio avec anciennes bornes 759 DO ij=ijb,ije760 761 if (q(ij,l,iq).gt.qperemin) then ! modif 13 nov 2020762 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)763 else764 Ratio(ij,l,iq2)=ratiomin765 endif751 DO ij=ijb,ije 752 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 753 if (q(ij,l,iq).gt.qperemin) then ! modif 13 nov 2020 754 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 755 else 756 Ratio(ij,l,iq2)=ratiomin 757 endif 766 758 enddo !DO ij=ijbm,ijem 767 enddo !DO l=1,llm 768 c$OMP END DO NOWAIT 769 enddo !do ifils=1,nqdesc(iq) 770 771 do ifils=1,nqfils(iq) 772 iq2=iqfils(ifils,iq) 773 call vly_loc(Ratio,pente_max,masse,qbyv,iq2) 774 enddo !do ifils=1,nqfils(iq) 775 endif !if (nqfils(iq).gt.0) then 759 enddo !DO l=1,llm 760 c$OMP END DO NOWAIT 761 enddo 762 763 do ifils=1,tracers(iq)%nqChilds 764 iq2=tracers(iq)%iqDescen(ifils) 765 call vly_loc(Ratio,pente_max,masse,qbyv,iq2) 766 enddo 776 767 ! end CRisi 777 768 … … 862 853 ! if (pole_sud) ije=ij_end 863 854 864 if (nqfils(iq).gt.0) then 865 do ifils=1,nqdesc(iq) 866 iq2=iqfils(ifils,iq) 855 do ifils=1,tracers(iq)%nqDescen 856 iq2=tracers(iq)%iqDescen(ifils) 867 857 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 868 858 DO l=1,llm 869 859 DO ij=ijb,ije 870 860 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 871 861 enddo 872 enddo 873 c$OMP END DO NOWAIT 874 enddo !do ifils=1,nqdesc(iq) 875 endif !if (nqfils(iq).gt.0) then 862 enddo 863 c$OMP END DO NOWAIT 864 enddo 876 865 877 866 … … 895 884 USE parallel_lmdz 896 885 USE vlz_mod 897 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils, ! CRisi &886 USE infotrac, ONLY : nqtot,tracers, ! CRisi & 898 887 & qperemin,masseqmin,ratiomin ! MVals et CRisi 899 888 … … 1159 1148 ! CRisi: appel récursif de l'advection sur les fils. 1160 1149 ! Il faut faire ça avant d'avoir mis à jour q et masse 1161 !write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq) 1162 if (nqfils(iq).gt.0) then 1163 do ifils=1,nqdesc(iq) 1164 !do ifils=1,nqfils(iq) ! modif C Risi 22 nov 2020 1165 iq2=iqfils(ifils,iq) 1166 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1167 DO l=1,llm 1150 !write(*,*) 'vlsplt 942: iq,nqChilds(iq)=',iq,tracers(iq)%nqChilds 1151 do ifils=1,tracers(iq)%nqDescen 1152 iq2=tracers(iq)%iqDescen(ifils) 1153 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1154 DO l=1,llm 1168 1155 DO ij=ijb,ije 1169 1156 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 1170 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)1171 if (q(ij,l,iq).gt.qperemin) then1172 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)1173 else1174 Ratio(ij,l,iq2)=ratiomin1175 endif1176 !wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai20151177 w(ij,l,iq2)=wq(ij,l,iq)1157 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin) 1158 if (q(ij,l,iq).gt.qperemin) then 1159 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 1160 else 1161 Ratio(ij,l,iq2)=ratiomin 1162 endif 1163 !wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015 1164 w(ij,l,iq2)=wq(ij,l,iq) 1178 1165 enddo 1179 1180 c$OMP END DO NOWAIT 1181 enddo !do ifils=1,nqdesc(iq)1166 enddo 1167 c$OMP END DO NOWAIT 1168 enddo 1182 1169 c$OMP BARRIER 1183 1170 1184 do ifils=1,nqfils(iq) 1185 iq2=iqfils(ifils,iq) 1186 call vlz_loc(Ratio,pente_max,masse,w,ijb_x,ije_x,iq2) 1187 enddo !do ifils=1,nqfils(iq) 1188 endif !if (nqfils(iq).gt.0) then 1171 do ifils=1,tracers(iq)%nqChilds 1172 iq2=tracers(iq)%iqDescen(ifils) 1173 call vlz_loc(Ratio,pente_max,masse,w,ijb_x,ije_x,iq2) 1174 enddo 1189 1175 ! end CRisi 1190 1176 … … 1207 1193 1208 1194 ! retablir les fils en rapport de melange par rapport a l'air: 1209 if (nqfils(iq).gt.0) then 1210 do ifils=1,nqdesc(iq) 1211 iq2=iqfils(ifils,iq) 1195 do ifils=1,tracers(iq)%nqDescen 1196 iq2=tracers(iq)%iqDescen(ifils) 1212 1197 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1213 1198 DO l=1,llm 1214 1199 DO ij=ijb,ije 1215 1200 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 1216 1201 enddo 1217 enddo 1218 c$OMP END DO NOWAIT 1219 enddo !do ifils=1,nqdesc(iq) 1220 endif !if (nqfils(iq).gt.0) then 1202 enddo 1203 c$OMP END DO NOWAIT 1204 enddo 1221 1205 1222 1206 RETURN -
LMDZ6/trunk/libf/dyn3dmem/vlspltgen_loc.F
r2603 r4050 2 2 ! $Header$ 3 3 ! 4 SUBROUTINE vlspltgen_loc( q, iadv,pente_max,masse,w,pbaru,pbarv,4 SUBROUTINE vlspltgen_loc( q,pente_max,masse,w,pbaru,pbarv, 5 5 & pdt, p,pk,teta ) 6 6 … … 28 28 USE VAMPIR 29 29 ! CRisi: on rajoute variables utiles d'infotrac 30 USE infotrac, ONLY : nqtot,nqperes,nqdesc,nqfils,iqfils, 31 & ok_iso_verif 30 USE infotrac, ONLY : nqtot,nqperes, tracers,ok_iso_verif 32 31 USE vlspltgen_mod 33 32 USE comconst_mod, ONLY: cpp … … 41 40 c Arguments: 42 41 c ---------- 43 INTEGER iadv(nqtot)44 42 REAL masse(ijb_u:ije_u,llm),pente_max 45 43 REAL pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm) … … 200 198 ! DO iq=1,nqtot 201 199 DO iq=1,nqperes ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air 202 !write(*,*) 'vlspltgen 192: iq,iadv=',iq, iadv(iq)200 !write(*,*) 'vlspltgen 192: iq,iadv=',iq,tracers(iq)%iadv 203 201 #ifdef DEBUG_IO 204 202 CALL WriteField_u('zq',zq(:,:,iq)) 205 203 CALL WriteField_u('zm',zm(:,:,iq)) 206 204 #endif 207 if(iadv(iq) == 0) then 208 209 cycle 210 211 else if (iadv(iq)==10) then 212 205 SELECT CASE(tracers(iq)%iadv) 206 CASE(0); CYCLE 207 CASE(10) 213 208 #ifdef _ADV_HALO 214 209 ! CRisi: on ajoute les nombres de fils et tableaux des fils … … 229 224 call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1) 230 225 ! CRisi 231 do ifils=1, nqdesc(iq)232 iq2= iqfils(ifils,iq)226 do ifils=1,tracers(iq)%nqDescen 227 iq2=tracers(iq)%iqDescen(ifils) 233 228 call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1) 234 229 call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1) … … 238 233 call VTe(VTHallo) 239 234 c$OMP END MASTER 240 else if (iadv(iq)==14) then 241 235 CASE(14) 242 236 #ifdef _ADV_HALO 243 237 call vlxqs_loc(zq,pente_max,zm,mu, … … 256 250 call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1) 257 251 call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1) 258 do ifils=1, nqdesc(iq)259 iq2= iqfils(ifils,iq)252 do ifils=1,tracers(iq)%nqDescen 253 iq2=tracers(iq)%iqDescen(ifils) 260 254 call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1) 261 255 call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1) … … 265 259 call VTe(VTHallo) 266 260 c$OMP END MASTER 267 else 268 261 CASE DEFAULT 269 262 stop 'vlspltgen_p : schema non parallelise' 270 263 271 endif264 END SELECT 272 265 273 266 enddo !DO iq=1,nqperes … … 298 291 !write(*,*) 'vlspltgen 279: iq=',iq 299 292 300 if(iadv(iq) == 0) then 301 302 cycle 303 304 else if (iadv(iq)==10) then 305 293 SELECT CASE(tracers(iq)%iadv) 294 CASE(0); CYCLE 295 CASE(10) 306 296 #ifdef _ADV_HALLO 307 297 call vlx_loc(zq,pente_max,zm,mu, 308 298 & ij_begin+2*iip1,ij_end-2*iip1,iq) 309 299 #endif 310 else if (iadv(iq)==14) then300 CASE(14) 311 301 #ifdef _ADV_HALLO 312 302 call vlxqs_loc(zq,pente_max,zm,mu, 313 303 & qsat,ij_begin+2*iip1,ij_end-2*iip1,iq) 314 304 #endif 315 else 316 305 CASE DEFAULT 317 306 stop 'vlspltgen_p : schema non parallelise' 318 307 319 endif308 END SELECT 320 309 321 310 enddo … … 355 344 #endif 356 345 357 if(iadv(iq) == 0) then 358 359 cycle 360 361 else if (iadv(iq)==10) then 362 363 call vly_loc(zq,pente_max,zm,mv,iq) 364 365 else if (iadv(iq)==14) then 366 367 call vlyqs_loc(zq,pente_max,zm,mv, 368 & qsat,iq) 369 370 else 371 346 SELECT CASE(tracers(iq)%iadv) 347 CASE(0); CYCLE 348 CASE(10); call vly_loc(zq,pente_max,zm,mv,iq) 349 CASE(14); call vlyqs_loc(zq,pente_max,zm,mv,qsat,iq) 350 CASE DEFAULT 372 351 stop 'vlspltgen_p : schema non parallelise' 373 374 endif 352 END SELECT 375 353 376 354 enddo … … 386 364 CALL WriteField_u('zm',zm(:,:,iq)) 387 365 #endif 388 if(iadv(iq) == 0) then 389 390 cycle 391 392 else if (iadv(iq)==10 .or. iadv(iq)==14 ) then 393 366 SELECT CASE(tracers(iq)%iadv) 367 CASE(0); CYCLE 368 CASE(10,14) 394 369 c$OMP BARRIER 395 370 #ifdef _ADV_HALLO … … 411 386 call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2) 412 387 ! CRisi 413 do ifils=1, nqdesc(iq)414 iq2= iqfils(ifils,iq)388 do ifils=1,tracers(iq)%nqDescen 389 iq2=tracers(iq)%iqDescen(ifils) 415 390 call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest2) 416 391 call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest2) … … 420 395 c$OMP END MASTER 421 396 c$OMP BARRIER 422 else 423 397 CASE DEFAULT 424 398 stop 'vlspltgen_p : schema non parallelise' 425 399 426 endif400 END SELECT 427 401 428 402 enddo … … 448 422 !write(*,*) 'vlspltgen 409: iq=',iq 449 423 450 if(iadv(iq) == 0) then 451 452 cycle 453 454 else if (iadv(iq)==10 .or. iadv(iq)==14 ) then 424 SELECT CASE(tracers(iq)%iadv) 425 CASE(0); CYCLE 426 CASE(10,14) 455 427 c$OMP BARRIER 456 428 … … 461 433 462 434 c$OMP BARRIER 463 else 464 435 CASE DEFAULT 465 436 stop 'vlspltgen_p : schema non parallelise' 466 467 endif 437 END SELECT 468 438 469 439 enddo … … 498 468 CALL WriteField_u('zm',zm(:,:,iq)) 499 469 #endif 500 if(iadv(iq) == 0) then 501 502 cycle 503 504 else if (iadv(iq)==10) then 505 506 call vly_loc(zq,pente_max,zm,mv,iq) 507 508 else if (iadv(iq)==14) then 509 510 call vlyqs_loc(zq,pente_max,zm,mv, 511 & qsat,iq) 512 513 else 514 515 stop 'vlspltgen_p : schema non parallelise' 516 517 endif 470 SELECT CASE(tracers(iq)%iadv) 471 CASE(0); CYCLE 472 CASE(10); call vly_loc(zq,pente_max,zm,mv, iq) 473 CASE(14); call vlyqs_loc(zq,pente_max,zm,mv,qsat,iq) 474 CASE DEFAULT; stop 'vlspltgen_p : schema non parallelise' 475 END SELECT 518 476 519 477 enddo !do iq=1,nqperes … … 529 487 CALL WriteField_u('zm',zm(:,:,iq)) 530 488 #endif 531 if(iadv(iq) == 0) then 532 533 cycle 534 535 else if (iadv(iq)==10) then 536 537 call vlx_loc(zq,pente_max,zm,mu, 489 SELECT CASE(tracers(iq)%iadv) 490 CASE(0); CYCLE 491 CASE(10); call vlx_loc(zq,pente_max,zm,mu, 538 492 & ij_begin,ij_end,iq) 539 540 else if (iadv(iq)==14) then 541 542 call vlxqs_loc(zq,pente_max,zm,mu, 493 CASE(14); call vlxqs_loc(zq,pente_max,zm,mu, 543 494 & qsat, ij_begin,ij_end,iq) 544 545 else 546 547 stop 'vlspltgen_p : schema non parallelise' 548 549 endif 495 CASE DEFAULT; stop 'vlspltgen_p : schema non parallelise' 496 END SELECT 550 497 551 498 enddo !do iq=1,nqperes -
LMDZ6/trunk/libf/dyn3dmem/vlspltqs_loc.F
r3800 r4050 12 12 c -------------------------------------------------------------------- 13 13 USE parallel_lmdz 14 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils, ! CRisi &14 USE infotrac, ONLY : nqtot,tracers, ! CRisi & 15 15 & qperemin,masseqmin,ratiomin ! MVals et CRisi 16 16 IMPLICIT NONE … … 337 337 ! CRisi: appel récursif de l'advection sur les fils. 338 338 ! Il faut faire ça avant d'avoir mis à jour q et masse 339 !write(*,*) 'vlspltqs 336: iq,ijb_x,nqfils(iq)=', 340 ! & iq,ijb_x,nqfils(iq) 341 342 if (nqfils(iq).gt.0) then 343 do ifils=1,nqdesc(iq) 344 iq2=iqfils(ifils,iq) 339 !write(*,*) 'vlspltqs 336: iq,ijb_x,iqDescen(iq)=', 340 ! & iq,ijb_x,tracers(iq)%iqDescen 341 342 do ifils=1,tracers(iq)%nqDescen 343 iq2=tracers(iq)%iqDescen(ifils) 345 344 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 346 345 DO l=1,llm 347 346 DO ij=ijb,ije 348 !MVals: veiller a ce qu'on n'ait pas de denominateur nul349 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)350 if (q(ij,l,iq).gt.qperemin) then ! modif 13 nov 2020351 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)352 else353 Ratio(ij,l,iq2)=ratiomin354 endif347 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 348 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin) 349 if (q(ij,l,iq).gt.qperemin) then ! modif 13 nov 2020 350 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 351 else 352 Ratio(ij,l,iq2)=ratiomin 353 endif 355 354 enddo 356 enddo 357 c$OMP END DO NOWAIT 358 enddo !do ifils=1,nqfils(iq) 359 do ifils=1,nqfils(iq) 360 iq2=iqfils(ifils,iq) 361 !write(*,*) 'vlxqs 349: on appelle vlx pour iq2=',iq2 362 call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2) 363 enddo !do ifils=1,nqfils(iq) 364 endif !if (nqfils(iq).gt.0) then 355 enddo 356 c$OMP END DO NOWAIT 357 enddo 358 do ifils=1,tracers(iq)%nqDescen 359 iq2=tracers(iq)%iqDescen(ifils) 360 !write(*,*) 'vlxqs 349: on appelle vlx pour iq2=',iq2 361 call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2) 362 enddo 365 363 ! end CRisi 366 364 … … 389 387 390 388 ! retablir les fils en rapport de melange par rapport a l'air: 391 if (nqfils(iq).gt.0) then 392 do ifils=1,nqdesc(iq) 393 iq2=iqfils(ifils,iq) 389 do ifils=1,tracers(iq)%nqDescen 390 iq2=tracers(iq)%iqDescen(ifils) 394 391 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 395 392 DO l=1,llm 396 393 DO ij=ijb+1,ije 397 394 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 398 395 enddo 399 396 DO ij=ijb+iip1-1,ije,iip1 400 397 q(ij-iim,l,iq2)=q(ij,l,iq2) 401 398 enddo ! DO ij=ijb+iip1-1,ije,iip1 402 enddo 403 c$OMP END DO NOWAIT 404 enddo !do ifils=1,nqdesc(iq) 405 endif !if (nqfils(iq).gt.0) then 399 enddo 400 c$OMP END DO NOWAIT 401 enddo 406 402 407 403 !write(*,*) 'vlspltqs 399: iq,ijb_x=',iq,ijb_x … … 426 422 c -------------------------------------------------------------------- 427 423 USE parallel_lmdz 428 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils, ! CRisi &424 USE infotrac, ONLY : nqtot,tracers, ! CRisi & 429 425 & qperemin,masseqmin,ratiomin ! MVals et CRisi 430 426 USE comconst_mod, ONLY: pi … … 733 729 ! CRisi: appel récursif de l'advection sur les fils. 734 730 ! Il faut faire ça avant d'avoir mis à jour q et masse 735 !write(*,*) 'vlyqs 689: iq, nqfils(iq)=',iq,nqfils(iq)731 !write(*,*) 'vlyqs 689: iq,iqDescen(iq)=',iq,tracers(iq)%iqDescen 736 732 737 733 ijb=ij_begin-2*iip1 … … 747 743 !write(lunout,*) 'ij_begin,ij_end=',ij_begin,ij_end 748 744 !write(lunout,*) 'pole_nord,pole_sud=',pole_nord,pole_sud 749 if (nqfils(iq).gt.0) then 750 do ifils=1,nqdesc(iq) 751 iq2=iqfils(ifils,iq) 745 do ifils=1,tracers(iq)%nqDescen 746 iq2=tracers(iq)%iqDescen(ifils) 752 747 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 753 748 DO l=1,llm 754 749 ! modif des bornes: CRisi 16 nov 2020 755 750 ! d'abord masse avec bornes corrigées 756 751 DO ij=ijbm,ijem 757 !MVals: veiller a ce qu'on n'ait pas de denominateur nul758 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)752 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 753 masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin) 759 754 enddo !DO ij=ijbm,ijem 760 755 761 756 ! ensuite Ratio avec anciennes bornes 762 757 DO ij=ijb,ije 763 !MVals: veiller a ce qu'on n'ait pas de denominateur nul764 !write(lunout,*) 'ij,l,q(ij,l,iq)=',ij,l,q(ij,l,iq)765 if (q(ij,l,iq).gt.qperemin) then ! modif 13 nov 2020766 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)767 else768 Ratio(ij,l,iq2)=ratiomin769 endif758 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 759 !write(lunout,*) 'ij,l,q(ij,l,iq)=',ij,l,q(ij,l,iq) 760 if (q(ij,l,iq).gt.qperemin) then ! modif 13 nov 2020 761 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 762 else 763 Ratio(ij,l,iq2)=ratiomin 764 endif 770 765 enddo !DO ij=ijbm,ijem 771 enddo !DO l=1,llm 772 c$OMP END DO NOWAIT 773 enddo !do ifils=1,nqdesc(iq) 774 do ifils=1,nqfils(iq) 775 iq2=iqfils(ifils,iq) 776 !write(lunout,*) 'vly: appel recursiv vly iq2=',iq2 777 call vly_loc(Ratio,pente_max,masse,qbyv,iq2) 778 enddo !do ifils=1,nqfils(iq) 779 endif !if (nqfils(iq).gt.0) then 766 enddo !DO l=1,llm 767 c$OMP END DO NOWAIT 768 enddo 769 do ifils=1,tracers(iq)%nqDescen 770 iq2=tracers(iq)%iqDescen(ifils) 771 !write(lunout,*) 'vly: appel recursiv vly iq2=',iq2 772 call vly_loc(Ratio,pente_max,masse,qbyv,iq2) 773 enddo 780 774 781 775 … … 856 850 ! if (pole_sud) ije=ij_end-iip1 857 851 858 if (nqfils(iq).gt.0) then 859 do ifils=1,nqdesc(iq) 860 iq2=iqfils(ifils,iq) 852 do ifils=1,tracers(iq)%nqDescen 853 iq2=tracers(iq)%iqDescen(ifils) 861 854 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 862 855 DO l=1,llm 863 856 DO ij=ijb,ije 864 857 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 865 858 enddo 866 enddo 867 c$OMP END DO NOWAIT 868 enddo !do ifils=1,nqdesc(iq) 869 endif !if (nqfils(iq).gt.0) then 859 enddo 860 c$OMP END DO NOWAIT 861 enddo 870 862 871 863 -
LMDZ6/trunk/libf/dyn3dmem/vlz_mod.F90
r2281 r4050 5 5 REAL,POINTER,SAVE :: dzqw(:,:) 6 6 REAL,POINTER,SAVE :: adzqw(:,:) 7 ! CRisi: pour les traceurs: 8 !REAL,POINTER,SAVE :: masseq(:,:,:) 7 ! CRisi: pour les traceurs: 9 8 REAL,POINTER,SAVE :: Ratio(:,:,:) 10 9 … … 25 24 CALL allocate_u(dzqw,llm,d) 26 25 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 26 IF(ANY(tracers(:)%nqDescen > 0)) CALL allocate_u(Ratio,llm,nqtot,d) 31 27 32 28 END SUBROUTINE vlz_allocate … … 44 40 CALL switch_u(dzqw,distrib_vanleer,dist) 45 41 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 42 IF(ANY(tracers(:)%nqDescen > 0)) CALL switch_u(Ratio,distrib_vanleer,dist) 51 43 52 44 END SUBROUTINE vlz_switch_vanleer -
LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0dyn_netcdf.F90
r4049 r4050 150 150 CALL regr_lat_time_coefoz 151 151 CALL press_coefoz 152 CALL regr_pr_o3(p3d, q3d(:,:,:,i ))153 q3d(:,:,:,i )=q3d(:,:,:,i)*48./ 29.!--- Mole->mass fraction152 CALL regr_pr_o3(p3d, q3d(:,:,:,iq)) 153 q3d(:,:,:,iq)=q3d(:,:,:,iq)*48./ 29. !--- Mole->mass fraction 154 154 END IF 155 155 #endif -
LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90
r4046 r4050 18 18 USE infotrac, ONLY: nqtot,nqo,nbtr,nqCO2,tracers,type_trac,& 19 19 niadv,conv_flg,pbl_flg,solsym,& 20 nqfils,nqdesc,nqdesc_tot,iqfils,iqpere,&21 20 ok_isotopes,ok_iso_verif,ok_isotrac,& 22 21 ok_init_iso,niso_possibles,tnat,& 23 alpha_ideal,use_iso,iqiso,iso_num,& 24 iso_indnum,zone_num,phase_num,& 22 alpha_ideal,use_iso,iqiso,iso_indnum,& 25 23 indnum_fn_num,index_trac,& 26 24 niso,ntraceurs_zone,ntraciso,nqtottr,itr_indice … … 148 146 CALL init_infotrac_phy(nqtot,nqo,nbtr,nqtottr,nqCO2,tracers,type_trac,& 149 147 niadv,conv_flg,pbl_flg,solsym,& 150 nqfils,nqdesc,nqdesc_tot,iqfils,iqpere,&151 148 ok_isotopes,ok_iso_verif,ok_isotrac,& 152 149 ok_init_iso,niso_possibles,tnat,& 153 alpha_ideal,use_iso,iqiso,iso_num,& 154 iso_indnum,zone_num,phase_num,& 150 alpha_ideal,use_iso,iqiso,iso_indnum,& 155 151 indnum_fn_num,index_trac,& 156 152 niso,ntraceurs_zone,ntraciso,itr_indice & -
LMDZ6/trunk/libf/phylmd/infotrac_phy.F90
r4046 r4050 51 51 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: niadv ! equivalent dyn / physique 52 52 !$OMP THREADPRIVATE(niadv) 53 54 ! CRisi: tableaux de fils55 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqfils56 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqdesc ! nombres de fils + nombre de tous les petits fils sur toutes les générations57 INTEGER, SAVE :: nqdesc_tot58 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqfils59 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iqpere60 !$OMP THREADPRIVATE(nqfils,nqdesc,nqdesc_tot,iqfils,iqpere)61 53 62 54 ! conv_flg(it)=0 : convection desactivated for tracer number it … … 84 76 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqiso ! donne indice iq en fn de (ixt,phase) 85 77 !$OMP THREADPRIVATE(iqiso) 86 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_num ! donne numéro iso entre 1 et niso_possibles en fn de nqtot87 !$OMP THREADPRIVATE(iso_num)88 78 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_indnum ! donne numéro iso entre 1 et niso effectif en fn de nqtot 89 79 !$OMP THREADPRIVATE(iso_indnum) 90 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: zone_num ! donne numéro de la zone de tracage en fn de nqtot91 !$OMP THREADPRIVATE(zone_num)92 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: phase_num ! donne numéro de la zone de tracage en fn de nqtot93 !$OMP THREADPRIVATE(phase_num)94 80 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 95 81 !$OMP THREADPRIVATE(indnum_fn_num) … … 106 92 SUBROUTINE init_infotrac_phy(nqtot_,nqo_,nbtr_,nqtottr_,nqCO2_,tracers_,type_trac_,& 107 93 niadv_,conv_flg_,pbl_flg_,solsym_,& 108 nqfils_,nqdesc_,nqdesc_tot_,iqfils_,iqpere_,&109 94 ok_isotopes_,ok_iso_verif_,ok_isotrac_,& 110 95 ok_init_iso_,niso_possibles_,tnat_,& 111 alpha_ideal_,use_iso_,iqiso_,iso_num_,& 112 iso_indnum_,zone_num_,phase_num_,& 96 alpha_ideal_,use_iso_,iqiso_,iso_indnum_,& 113 97 indnum_fn_num_,index_trac_,& 114 98 niso_,ntraceurs_zone_,ntraciso_,itr_indice_& … … 143 127 CHARACTER(len=*),INTENT(IN) :: solsym_(nbtr_) 144 128 ! Isotopes: 145 INTEGER,INTENT(IN) :: nqfils_(nqtot_)146 INTEGER,INTENT(IN) :: nqdesc_(nqtot_)147 INTEGER,INTENT(IN) :: nqdesc_tot_148 INTEGER,INTENT(IN) :: iqfils_(nqtot_,nqtot_)149 INTEGER,INTENT(IN) :: iqpere_(nqtot_)150 129 LOGICAL,INTENT(IN) :: ok_isotopes_ 151 130 LOGICAL,INTENT(IN) :: ok_iso_verif_ … … 157 136 LOGICAL,INTENT(IN) :: use_iso_(niso_possibles_) 158 137 INTEGER,INTENT(IN) :: iqiso_(ntraciso_,nqo_) 159 INTEGER,INTENT(IN) :: iso_num_(nqtot_)160 138 INTEGER,INTENT(IN) :: iso_indnum_(nqtot_) 161 INTEGER,INTENT(IN) :: zone_num_(nqtot_)162 INTEGER,INTENT(IN) :: phase_num_(nqtot_)163 139 INTEGER,INTENT(IN) :: indnum_fn_num_(niso_possibles_) 164 140 INTEGER,INTENT(IN) :: index_trac_(ntraceurs_zone_,niso_) … … 169 145 170 146 CHARACTER(LEN=30) :: modname="init_infotrac_phy" 147 INTEGER :: iq 171 148 172 149 nqtot=nqtot_ … … 176 153 nqtottr=nqtottr_ 177 154 ALLOCATE(tracers(nqtot)); tracers(:) = tracers_(:) 155 tracers(:)%isAdvected = tracers(:)%iadv > 0 156 ! tracers(:)%isH2Ofamily = delPhase(tracers(:)%gen0Name) == 'H2O' 157 tracers(:)%isH2Ofamily = [(tracers(iq)%gen0Name(1:3) == 'H2O', iq=1, nqtot)] 178 158 #ifdef CPP_StratAer 179 159 nbtr_bin=nbtr_bin_ … … 216 196 217 197 IF (ok_isotopes) THEN 218 ALLOCATE(nqfils(nqtot))219 nqfils(:)=nqfils_(:)220 ALLOCATE(nqdesc(nqtot))221 nqdesc(:)=nqdesc_(:)222 nqdesc_tot=nqdesc_tot_223 ALLOCATE(iqfils(nqtot,nqtot))224 iqfils(:,:)=iqfils_(:,:)225 ALLOCATE(iqpere(nqtot))226 iqpere(:)=iqpere_(:)227 228 198 tnat(:)=tnat_(:) 229 199 alpha_ideal(:)=alpha_ideal_(:) … … 232 202 ALLOCATE(iqiso(ntraciso,nqo)) 233 203 iqiso(:,:)=iqiso_(:,:) 234 ALLOCATE(iso_num(nqtot))235 iso_num(:)=iso_num_(:)236 204 ALLOCATE(iso_indnum(nqtot)) 237 205 iso_indnum(:)=iso_indnum_(:) 238 ALLOCATE(zone_num(nqtot))239 zone_num(:)=zone_num_(:)240 ALLOCATE(phase_num(nqtot))241 phase_num(:)=phase_num_(:)242 206 243 207 indnum_fn_num(:)=indnum_fn_num_(:) -
LMDZ6/trunk/libf/phylmd/tracco2i_mod.F90
r3857 r4050 30 30 31 31 USE dimphy 32 USE infotrac_phy 32 USE infotrac_phy, ONLY: nbtr 33 33 USE geometry_mod, ONLY: cell_area 34 34 USE carbon_cycle_mod, ONLY: id_CO2, nbcf_in, fields_in, cfname_in … … 336 336 337 337 USE dimphy 338 USE infotrac_phy338 ! USE infotrac_phy 339 339 USE geometry_mod, ONLY : cell_area 340 340 USE mod_grid_phy_lmdz -
LMDZ6/trunk/libf/phylmd/traclmdz_mod.F90
r4046 r4050 67 67 68 68 USE dimphy 69 USE infotrac_phy 69 USE infotrac_phy, ONLY: nbtr, tracers, niadv, solsym 70 70 71 71 ! Input argument … … 89 89 ! Initialization of the tracers should be done here only for those not found in the restart file. 90 90 USE dimphy 91 USE infotrac_phy 91 USE infotrac_phy, ONLY: nbtr, nqo, tracers, pbl_flg, conv_flg, niadv 92 92 USE regr_pr_comb_coefoz_m, ONLY: alloc_coefoz 93 93 USE press_coefoz_m, ONLY: press_coefoz … … 176 176 !! iiq=niadv(it+2) ! jyg 177 177 iiq=niadv(it+nqo) ! jyg 178 print*,'###'//TRIM(tracers(iiq)%name)//'###'179 print*,'###'//TRIM(strLower(tracers(iiq)%name))//'###'180 178 SELECT CASE(strLower(tracers(iiq)%name)) 181 179 CASE("rn"); id_rn = it ! radon … … 311 309 312 310 USE dimphy 313 USE infotrac_phy 311 USE infotrac_phy, ONLY: nbtr, pbl_flg, solsym 314 312 USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz 315 313 USE o3_chem_m, ONLY: o3_chem … … 586 584 ! variable trs is written to restart file (restartphy.nc) 587 585 USE dimphy 588 USE infotrac_phy 586 USE infotrac_phy, ONLY: nbtr 589 587 590 588 REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out -
LMDZ6/trunk/libf/phylmd/tracreprobus_mod.F90
r3666 r4050 12 12 13 13 USE dimphy 14 USE infotrac_phy 14 USE infotrac_phy, ONLY: nbtr, solsym 15 15 #ifdef REPROBUS 16 16 USE CHEM_REP, ONLY : pdt_rep, & ! pas de temps reprobus -
LMDZ6/trunk/libf/phylmdiso/cv30_routines.F90
r4004 r4050 449 449 450 450 #ifdef ISOVERIF 451 write(*,*) 'cv30_routine undilute 1 413: entr ée'451 write(*,*) 'cv30_routine undilute 1 413: entree' 452 452 #endif 453 453 … … 602 602 zfice(i) = MIN(MAX(zfice(i),0.0),1.0) 603 603 enddo 604 ! calcul de la composition du condensat glac éet liquide604 ! calcul de la composition du condensat glace et liquide 605 605 606 606 do i=1,len … … 647 647 648 648 #ifdef ISOVERIF 649 write(*,*) 'cv30_routine undilute 1 598: apr ès condiso'649 write(*,*) 'cv30_routine undilute 1 598: apres condiso' 650 650 651 651 if (iso_eau.gt.0) then … … 1012 1012 else 1013 1013 q(i,k)=0.0 1014 clw(i,k)=0.0 ! mise en commentaire le 5 avril pour v érif1014 clw(i,k)=0.0 ! mise en commentaire le 5 avril pour verif 1015 1015 ! convergence 1016 1016 endif !f (negation(essai_convergence)) then … … 1908 1908 real xtrti(ntraciso,nloc) 1909 1909 real xtres(ntraciso) 1910 ! on ajoute la dimension nloc à xtrti pour vérifs dans les tags: 5 fev1910 ! on ajoute la dimension nloc a xtrti pour verifs dans les tags: 5 fev 1911 1911 ! 2010 1912 1912 real zfice(nloc),zxtliq(ntraciso,nloc),zxtice(ntraciso,nloc) … … 1923 1923 #ifdef ISO 1924 1924 #ifdef ISOVERIF 1925 write(*,*) 'cv30_routines 1820: entr ée dans cv3_mixing'1925 write(*,*) 'cv30_routines 1820: entree dans cv3_mixing' 1926 1926 if (iso_eau.gt.0) then 1927 1927 call iso_verif_egalite_vect2D( & … … 1965 1965 xtelij(ixt,i,k,j)=0.0 1966 1966 enddo !do ixt =1,niso 1967 ! on initialise mieux que ça qent et elij, même si au final les1968 ! valeurs en nd=nl+1 ne sont pas utilis ées1967 ! on initialise mieux que ca qent et elij, meme si au final les 1968 ! valeurs en nd=nl+1 ne sont pas utilisees 1969 1969 qent(i,k,j)=rr(i,j) 1970 1970 elij(i,k,j)=0.0 … … 2121 2121 ! : 'tcond(il),rs(il,j)=', 2122 2122 ! : il,i,j,xtent(:,il,i,j),tcond(il),rs(il,j) 2123 ! colorier la vapeur r ésiduelle selon température de2124 ! condensation, et le condensat en un tag sp écifique2123 ! colorier la vapeur residuelle selon temperature de 2124 ! condensation, et le condensat en un tag spEcifique 2125 2125 if ((elij(il,i,j).gt.0.0).and.(qent(il,i,j).gt.0.0)) then 2126 2126 if (option_traceurs.eq.17) then … … 2241 2241 #ifdef ISOTRAC 2242 2242 if (option_tmin.ge.1) then 2243 ! colorier la vapeur r ésiduelle selon température de2244 ! condensation, et le condensat en un tag sp écifique2243 ! colorier la vapeur residuelle selon temperature de 2244 ! condensation, et le condensat en un tag specifique 2245 2245 ! write(*,*) 'cv3 tmp 2095 il,i,j,xtent(:,il,i,j)=', 2246 2246 ! : il,i,j,xtent(:,il,i,j) … … 2475 2475 #ifdef ISOTRAC 2476 2476 if (option_tmin.ge.1) then 2477 ! colorier la vapeur r ésiduelle selon température de2478 ! condensation, et le condensat en un tag sp écifique2477 ! colorier la vapeur residuelle selon temperature de 2478 ! condensation, et le condensat en un tag specifique 2479 2479 ! write(*,*) 'cv3 tmp 2314 il,i,j,xtent(:,il,i,j)=', 2480 2480 ! : il,i,j,xtent(:,il,i,j) … … 2579 2579 #ifdef ISO 2580 2580 #ifdef ISOTRAC 2581 ! seulement àla fin on taggue le condensat2581 ! seulement a la fin on taggue le condensat 2582 2582 if (option_cond.ge.1) then 2583 2583 do im = 1, nd 2584 2584 do jm = 1, nd 2585 2585 do il = 1, ncum 2586 ! colorier le condensat en un tag sp écifique2586 ! colorier le condensat en un tag specifique 2587 2587 do ixt=niso+1,ntraciso 2588 2588 if (index_zone(ixt).eq.izone_cond) then … … 2603 2603 do im = 1, nd 2604 2604 do il = 1, ncum 2605 ! colorier le condensat en un tag sp écifique2605 ! colorier le condensat en un tag specifique 2606 2606 do ixt=niso+1,ntraciso 2607 2607 if (index_zone(ixt).eq.izone_cond) then … … 2739 2739 ! ------------------------------------------------------ 2740 2740 !#ifdef ISOVERIF 2741 ! write(*,*) 'cv30_routines 2382: entr ée dans cv3_unsat'2741 ! write(*,*) 'cv30_routines 2382: entree dans cv3_unsat' 2742 2742 !#endif 2743 2743 … … 2747 2747 mp(:, :) = 0. 2748 2748 #ifdef ISO 2749 ! initialisation plus compl ète de water et rp2749 ! initialisation plus complete de water et rp 2750 2750 water(:,:)=0.0 2751 2751 xtwater(:,:,:)=0.0 … … 2936 2936 call iso_verif_traceur(xtwdtrain(1,il),'cv30_routine 2540') 2937 2937 if (option_cond.ge.1) then 2938 ! on v érifie que tout le détrainement est taggécondensat2938 ! on verifie que tout le detrainement est tagge condensat 2939 2939 if (iso_verif_positif_nostop( & 2940 2940 & xtwdtrain(index_trac(izone_cond,iso_eau),il) & … … 3032 3032 3033 3033 #ifdef ISO 3034 ! ajout cam: éviter les evaporations ou eaux négatives3035 ! water(il,i)=max(0.0,water(il,i)) ! ceci est toujours v érifié3034 ! ajout cam: eviter les evaporations ou eaux negatives 3035 ! water(il,i)=max(0.0,water(il,i)) ! ceci est toujours verifie 3036 3036 #ifdef ISOVERIF 3037 3037 call iso_verif_positif(water(il,i),'cv30_unsat 2376') … … 3189 3189 #ifdef ISO 3190 3190 #ifdef ISOVERIF 3191 ! verif des inputs àappel stewart3191 ! verif des inputs a appel stewart 3192 3192 ! write(*,*) 'cv30_routines 2842 tmp: appel de appel_stewart' 3193 3193 do il=1,ncum … … 3208 3208 enddo 3209 3209 #endif 3210 ! appel de appel_stewart_vectoris é3210 ! appel de appel_stewart_vectorise 3211 3211 call appel_stewart_vectall(lwork,ncum, & 3212 3212 & ph,t,evap,xtwdtrain, & … … 3268 3268 #endif 3269 3269 3270 ! équivalent isotopique de rp(il,i)=amin1(rp(il,i),rs(il,i))3270 ! equivalent isotopique de rp(il,i)=amin1(rp(il,i),rs(il,i)) 3271 3271 do il=1,ncum 3272 3272 if (i.lt.inb(il) .and. lwork(il)) then … … 3463 3463 real xtbx(ntraciso), xtawat(ntraciso) 3464 3464 ! cam debug 3465 ! pour l'homog énéisation sous le nuage:3465 ! pour l'homogeneisation sous le nuage: 3466 3466 real frsum(nloc), bxtsum(ntraciso,nloc), fxtsum(ntraciso,nloc) 3467 ! correction dans calcul tendance li ée àAm:3467 ! correction dans calcul tendance liee a Am: 3468 3468 real dq_tmp,k_tmp,dx_tmp,R_tmp,dqreste_tmp,dxreste_tmp,kad_tmp 3469 3469 logical correction_excess_aberrant 3470 3470 parameter (correction_excess_aberrant=.false.) 3471 ! correction qui permettait d' éviter deltas et dexcess aberrants. Mais3471 ! correction qui permettait d'eviter deltas et dexcess aberrants. Mais 3472 3472 ! pb: ne conserve pas la masse d'isotopes! 3473 3473 #ifdef DIAGISO 3474 ! diagnostiques juste: tendance des diff érents processus3474 ! diagnostiques juste: tendance des differents processus 3475 3475 real fxt_detrainement(ntraciso,nloc,nd) 3476 3476 real fxt_fluxmasse(ntraciso,nloc,nd) … … 3517 3517 #ifdef ISO 3518 3518 ! cam debug 3519 ! write(*,*) 'cv30_routines 3082: entr ée dans cv3_yield'3519 ! write(*,*) 'cv30_routines 3082: entree dans cv3_yield' 3520 3520 ! en cam debug 3521 3521 do ixt = 1, ntraciso … … 3749 3749 do ixt = 1, ntraciso 3750 3750 ! fxt_fluxmasse(ixt,il,1)=fxt_fluxmasse(ixt,il,1) & 3751 ! & +0.01*grav*am(il)*(xt(ixt,il,2)-xt(ixt,il,1))*work(il) ! d éplacé3752 ! plus haut car il existe diff érents cas3751 ! & +0.01*grav*am(il)*(xt(ixt,il,2)-xt(ixt,il,1))*work(il) ! deplace 3752 ! plus haut car il existe differents cas 3753 3753 fxt_ddft(ixt,il,1)=fxt_ddft(ixt,il,1) & 3754 3754 & +0.01*grav*mp(il,2)*(xtp(ixt,il,2)-xt(ixt,il,1))*work(il) … … 3759 3759 3760 3760 3761 ! pour l'ajout de la tendance li ée au flux de masse Am, il faut être3761 ! pour l'ajout de la tendance liee au flux de masse Am, il faut etre 3762 3762 ! prudent. 3763 3763 ! On a dq1=k*(q2-q1) avec k=dt*0.01*grav*am(il)*work(il) 3764 ! Pour les isotopes, la formule utilis ée depuis 2006 et qui avait toujours marchéest:3764 ! Pour les isotopes, la formule utilisee depuis 2006 et qui avait toujours marche est: 3765 3765 ! dx1=k*(x2-x1) 3766 ! Mais on plante dans un cas pathologique en d écembre 2017 lors du test3767 ! d'un cas d'Anne Cozic: les isotopes deviennent n égatifs.3766 ! Mais on plante dans un cas pathologique en decembre 2017 lors du test 3767 ! d'un cas d'Anne Cozic: les isotopes deviennent negatifs. 3768 3768 ! C'est un cas pas physique: on perd 99% de la masse de vapeur d'eau! 3769 3769 ! q2=1.01e-3 et q1=1.25e-3 kg/kg 3770 ! et dq=-1.24e-3: comment est-ce possible qu'un flux venant d'un air à3771 ! q2= 1.01e-3 ass èche q1 jusqu'à0.01e-3kg/kg!3772 ! Pour les isotopes, ça donne des x1+dx négatifs.3773 ! Ce n'est pas physique mais il faut quand m ême s'adapter.3774 ! Pour cela, on consid ère que d'abord on fait rentrer le flux de masse3770 ! et dq=-1.24e-3: comment est-ce possible qu'un flux venant d'un air a 3771 ! q2= 1.01e-3 asseche q1 jusqu'a 0.01e-3kg/kg! 3772 ! Pour les isotopes, ca donne des x1+dx negatifs. 3773 ! Ce n'est pas physique mais il faut quand meme s'adapter. 3774 ! Pour cela, on considere que d'abord on fait rentrer le flux de masse 3775 3775 ! descendant, et ensuite seulement on fait sortir le flux de masse 3776 3776 ! sortant. … … 3778 3778 ! isotopique de la vapeur d'eau q1. 3779 3779 ! A la fin, on a R=(x1+dx)/(q1+dq)=(x1+k*x2)/(q1+k*q2) 3780 ! On v érifie que quand k est petit, on tend vers la formulation3780 ! On verifie que quand k est petit, on tend vers la formulation 3781 3781 ! habituelle. 3782 ! Comme on est habitu és àla formulation habituelle, qu'elle a fait ses3783 ! preuves, on la garde sauf dans le cas o ù dq/q<-0.9 oùon utilise la3782 ! Comme on est habitues a la formulation habituelle, qu'elle a fait ses 3783 ! preuves, on la garde sauf dans le cas ou dq/q<-0.9 ou on utilise la 3784 3784 ! nouvelle formulation. 3785 3785 ! rappel: dq_tmp=0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)*delt 3786 ! M ême avec cette nouvelle foirmulation, on a encore des isotopes3787 ! n égatifs, cette fois àcause des ddfts3788 ! On consid ère donc les tendances et série et non en parallèle quand on3786 ! Meme avec cette nouvelle foirmulation, on a encore des isotopes 3787 ! negatifs, cette fois a cause des ddfts 3788 ! On considere donc les tendances et serie et non en parallele quand on 3789 3789 ! calcule R_tmp. 3790 3790 dq_tmp=0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)*delt ! utile ci-dessous 3791 3791 if ((dq_tmp/rr(il,1).lt.-0.9).and.correction_excess_aberrant) then 3792 ! nouvelle formulation o ùon fait d'abord entrer k*q2 et ensuite3792 ! nouvelle formulation ou on fait d'abord entrer k*q2 et ensuite 3793 3793 ! seulement on fait sortir k*q1 sans changement de composition 3794 3794 ! isotopique … … 3828 3828 enddo ! do ixt = 1, ntraciso 3829 3829 else !if (dq_tmp/rr(il,1).lt.-0.9) then 3830 ! formulation habituelle qui avait toujours march é de 2006 à3831 ! d écembre 2017.3830 ! formulation habituelle qui avait toujours marche de 2006 a 3831 ! decembre 2017. 3832 3832 do ixt = 1, ntraciso 3833 3833 fxt(ixt,il,1)=fxt(ixt,il,1) & … … 4232 4232 ! ad. 4233 4233 #endif 4234 ! ici, on s épare 2 cas, pour éviter le cas pathologique décrit plus haut4235 ! pour la tendance li ée à Am en i=1, qui peut conduire àdes isotopes4236 ! n égatifs dans les cas oùles flux de masse soustrait plus de 90% de la4237 ! vapeur de la couche. Voir plus haut le d étail des équations.4238 ! La diff érence ici est qu'on considère les flux de masse amp1 et ad en4239 ! m ême temps.4234 ! ici, on separe 2 cas, pour eviter le cas pathologique decrit plus haut 4235 ! pour la tendance liee a Am en i=1, qui peut conduire a des isotopes 4236 ! negatifs dans les cas ou les flux de masse soustrait plus de 90% de la 4237 ! vapeur de la couche. Voir plus haut le detail des equations. 4238 ! La difference ici est qu'on considere les flux de masse amp1 et ad en 4239 ! meme temps. 4240 4240 dq_tmp= 0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) & 4241 4241 & -ad(il)*(rr(il,i)-rr(il,i-1)))*delt 4242 ! c'est équivalent àdqi= kamp1*qip1+kad*qim1-(kamp1+kad)*qi4242 ! c'est equivalent a dqi= kamp1*qip1+kad*qim1-(kamp1+kad)*qi 4243 4243 if ((dq_tmp/rr(il,i).lt.-0.9).and.correction_excess_aberrant) then 4244 4244 ! nouvelle formulation … … 4430 4430 ! on change le traitement de cette ligne le 8 mai 2009: 4431 4431 ! avant, on avait: xtawat=xtelij(il,k,i)-(1.-xtep(il,i))*xtclw(il,i) 4432 ! c'est àdire que Rawat=Relij+(1-ep)*clw/awat*(Relij-Rclw)4433 ! si Relij!=Rclw, alors un fractionnement isotopique non physique était4432 ! c'est a dire que Rawat=Relij+(1-ep)*clw/awat*(Relij-Rclw) 4433 ! si Relij!=Rclw, alors un fractionnement isotopique non physique etait 4434 4434 ! introduit. 4435 ! En fait, awat repr ésente le surplus de condensat dans le mélange par4436 ! rapport àcelui restant dans la colonne adiabatique4437 ! ce surplus à la même compo que le elij, sans fractionnement.4438 ! d'o ùle nouveau traitement ci-dessous.4435 ! En fait, awat represente le surplus de condensat dans le melange par 4436 ! rapport a celui restant dans la colonne adiabatique 4437 ! ce surplus a la meme compo que le elij, sans fractionnement. 4438 ! d'ou le nouveau traitement ci-dessous. 4439 4439 if (elij(il,k,i).gt.0.0) then 4440 4440 do ixt = 1, ntraciso 4441 4441 xtawat(ixt)=awat*(xtelij(ixt,il,k,i)/elij(il,k,i)) 4442 ! xtawat(ixt)=amax1(xtawat(ixt),0.0) ! pas n écessaire4442 ! xtawat(ixt)=amax1(xtawat(ixt),0.0) ! pas necessaire 4443 4443 enddo 4444 4444 else !if (elij(il,k,i).gt.0.0) then 4445 4445 ! normalement, si elij(il,k,i)<=0, alors awat=0 4446 ! on le v érifie. Si c'est vrai -> xtawat=0 aussi4446 ! on le verifie. Si c'est vrai -> xtawat=0 aussi 4447 4447 #ifdef ISOVERIF 4448 4448 call iso_verif_egalite(awat,0.0,'cv30_yield 3779') … … 4942 4942 & 'cv30_yield 5029,O18, evap') 4943 4943 if ((il.eq.1636).and.(i.eq.9)) then 4944 write(*,*) 'cv30_yield 5057: ici, on v érifie deltaD_nobx'4944 write(*,*) 'cv30_yield 5057: ici, on verifie deltaD_nobx' 4945 4945 write(*,*) 'il,i=',il,i 4946 4946 write(*,*) 'fr(il,i),bx,fr(il,i)-bx=',fr(il,i),bx,fr(il,i)-bx … … 4973 4973 else ! taggage des ddfts: 4974 4974 ! la formule pour fq_ddft suppose que le ddft est en RP. Ce n'est pas le 4975 ! cas pour le water tagging puisqu'il y a conversion des mol écules4976 ! blances entrain ées en molécule rouges.4975 ! cas pour le water tagging puisqu'il y a conversion des molecules 4976 ! blances entrainees en molecule rouges. 4977 4977 ! Il faut donc prendre en compte ce taux de conversion quand 4978 4978 ! entrainement d'env vers ddft … … 4983 4983 ! : -conversion(iiso) 4984 4984 4985 ! Pb: quand on discretise, dqp/dt n'est pas v érifée numériquement.4986 ! on se retrouve donc avec des d Ye/dt diff érents de 0 même si ye=0 ( on4987 ! note X les mol écules poubelles et Y les molécules ddfts).4985 ! Pb: quand on discretise, dqp/dt n'est pas verifee numeriquement. 4986 ! on se retrouve donc avec des d Ye/dt differents de 0 meme si ye=0 ( on 4987 ! note X les molecules poubelles et Y les molecules ddfts). 4988 4988 4989 4989 ! Solution alternative: Dans le cas entrainant, Ye ne varie que par 4990 4990 ! ascendance compensatoire des ddfts et par perte de Ye vers le ddft. On 4991 ! calcule donc ce terme directement avec sch éma amont:4992 4993 ! ajout d éjà de l'évap4991 ! calcule donc ce terme directement avec schema amont: 4992 4993 ! ajout deja de l'evap 4994 4994 do ixt = 1+niso,ntraciso 4995 4995 fxt(ixt,il,i)=fxt(ixt,il,i) & … … 5069 5069 #endif 5070 5070 else !if (abs(dXe).gt.ridicule) then 5071 ! dans ce cas, fxtXe doit être faible5071 ! dans ce cas, fxtXe doit etre faible 5072 5072 5073 5073 #ifdef ISOVERIF … … 5085 5085 fxt(ixt,il,i)=fxt(ixt,il,i)+fxtXe(iiso) 5086 5086 else !if (izone.eq.izone_poubelle) then 5087 ! pas de tendance pour ce tag l à5087 ! pas de tendance pour ce tag la 5088 5088 endif !if (izone.eq.izone_poubelle) then 5089 5089 endif !if ((izone.ne.izone_revap).and. … … 5099 5099 5100 5100 else !if (mp(il,i).gt.mp(il,i+1)) then 5101 ! cas d étrainant: pas de problèmes5101 ! cas detrainant: pas de problemes 5102 5102 do ixt=1+niso,ntraciso 5103 5103 fxt(ixt,il,i)=fxt(ixt,il,i) & … … 5389 5389 DO il = 1, ncum 5390 5390 5391 ! attention, on corrige un probl ème C Risi5391 ! attention, on corrige un probleme C Risi 5392 5392 IF (cvflag_grav) then 5393 5393 … … 5722 5722 ! write(*,*) 'cv30_routine 3990: fin des il pour i=',i 5723 5723 enddo !do i=1,nl 5724 ! write(*,*) 'cv30_routine 3990: fin des v érifs sur homogen'5724 ! write(*,*) 'cv30_routine 3990: fin des verifs sur homogen' 5725 5725 #endif 5726 5726 … … 6027 6027 6028 6028 ! fraction deau condensee dans les melanges convertie en precip : epm 6029 ! et eau condens ée précipitée dans masse d'air saturé: l_m*dM_m/dzdz.dzdz6029 ! et eau condensee precipitee dans masse d'air sature : l_m*dM_m/dzdz.dzdz 6030 6030 DO j = 1, nam1 6031 6031 DO k = 1, j - 1 … … 6226 6226 6227 6227 #ifdef ISOVERIF 6228 write(*,*) 'cv30_routines 4293: entr ée dans cv3_uncompress'6228 write(*,*) 'cv30_routines 4293: entree dans cv3_uncompress' 6229 6229 #endif 6230 6230 DO i = 1, ncum … … 6346 6346 6347 6347 ! On fait varier epmax en fn de la cape 6348 ! Il faut donc recalculer ep, et hp qui a d éjà été calculéet6349 ! qui en d épend6350 ! Toutes les autres variables fn de ep sont calcul ées plus bas.6348 ! Il faut donc recalculer ep, et hp qui a deja ete calcule et 6349 ! qui en depend 6350 ! Toutes les autres variables fn de ep sont calculees plus bas. 6351 6351 6352 6352 #include "cvthermo.h" -
LMDZ6/trunk/libf/phylmdiso/cv3_routines.F90
r4033 r4050 403 403 enddo !do i=1,len 404 404 #endif 405 ! initialiser quelques variables oubli ées405 ! initialiser quelques variables oubliees 406 406 do i=1,len 407 407 plcllo(i)=0.0 … … 900 900 zfice(i) = MIN(MAX(zfice(i),0.0),1.0) 901 901 enddo 902 ! calcul de la composition du condensat glac éet liquide902 ! calcul de la composition du condensat glace et liquide 903 903 904 904 do i=1,len … … 959 959 960 960 #ifdef ISOVERIF 961 write(*,*) 'cv3_routine undilute 1 598: apr ès condiso'961 write(*,*) 'cv3_routine undilute 1 598: apres condiso' 962 962 963 963 if (iso_eau.gt.0) then … … 1435 1435 1436 1436 !JAM-------------------------------------------------------------------- 1437 ! Calcul de la quantit éd'eau sous forme de glace1437 ! Calcul de la quantite d'eau sous forme de glace 1438 1438 ! -------------------------------------------------------------------- 1439 1439 INTEGER nl, len … … 2856 2856 real xtrti(ntraciso,nloc) 2857 2857 real xtres(ntraciso) 2858 ! on ajoute la dimension nloc à xtrti pour vérifs dans les tags: 5 fev2858 ! on ajoute la dimension nloc a xtrti pour verifs dans les tags: 5 fev 2859 2859 ! 2010 2860 2860 real zfice(nloc),zxtliq(ntraciso,nloc),zxtice(ntraciso,nloc) … … 2873 2873 #ifdef ISO 2874 2874 #ifdef ISOVERIF 2875 ! write(*,*) 'cv3_routines 1820: entr ée dans cv3_mixing'2875 ! write(*,*) 'cv3_routines 1820: entree dans cv3_mixing' 2876 2876 do i=minorig+1,nl 2877 2877 do il=1,ncum … … 3083 3083 ! : 'tcond(il),rs(il,j)=', 3084 3084 ! : il,i,j,xtent(:,il,i,j),tcond(il),rs(il,j) 3085 ! colorier la vapeur r ésiduelle selon température de3086 ! condensation, et le condensat en un tag sp écifique3085 ! colorier la vapeur residuelle selon temperature de 3086 ! condensation, et le condensat en un tag specifique 3087 3087 if ((elij(il,i,j).gt.0.0).and.(qent(il,i,j).gt.0.0)) then 3088 3088 if (option_traceurs.eq.17) then … … 3194 3194 #ifdef ISOTRAC 3195 3195 if (option_tmin.ge.1) then 3196 ! colorier la vapeur r ésiduelle selon température de3197 ! condensation, et le condensat en un tag sp écifique3196 ! colorier la vapeur residuelle selon temperature de 3197 ! condensation, et le condensat en un tag specifique 3198 3198 ! write(*,*) 'cv3 tmp 2095 il,i,j,xtent(:,il,i,j)=', 3199 3199 ! : il,i,j,xtent(:,il,i,j) … … 3431 3431 #ifdef ISOTRAC 3432 3432 if (option_tmin.ge.1) then 3433 ! colorier la vapeur r ésiduelle selon température de3434 ! condensation, et le condensat en un tag sp écifique3433 ! colorier la vapeur residuelle selon temperature de 3434 ! condensation, et le condensat en un tag specifique 3435 3435 ! write(*,*) 'cv3 tmp 2314 il,i,j,xtent(:,il,i,j)=', 3436 3436 ! : il,i,j,xtent(:,il,i,j) … … 3543 3543 #ifdef ISO 3544 3544 #ifdef ISOTRAC 3545 ! seulement àla fin on taggue le condensat3545 ! seulement a la fin on taggue le condensat 3546 3546 if (option_cond.ge.1) then 3547 3547 do im = 1, nd 3548 3548 do jm = 1, nd 3549 3549 do il = 1, ncum 3550 ! colorier le condensat en un tag sp écifique3550 ! colorier le condensat en un tag specifique 3551 3551 do ixt=niso+1,ntraciso 3552 3552 if (index_zone(ixt).eq.izone_cond) then … … 3567 3567 do im = 1, nd 3568 3568 do il = 1, ncum 3569 ! colorier le condensat en un tag sp écifique3569 ! colorier le condensat en un tag specifique 3570 3570 do ixt=niso+1,ntraciso 3571 3571 if (index_zone(ixt).eq.izone_cond) then … … 3991 3991 call iso_verif_traceur(xtwdtrain(1,il),'cv3_routine 2540') 3992 3992 if (option_cond.ge.1) then 3993 ! on v érifie que tout le détrainement est taggécondensat3993 ! on verifie que tout le detrainement est tagge condensat 3994 3994 if (iso_verif_positif_nostop( & 3995 3995 & xtwdtrain(index_trac(izone_cond,iso_eau),il) & … … 4156 4156 !!---end jyg--- 4157 4157 4158 ! --------retour àla formulation originale d'Emanuel.4158 ! --------retour a la formulation originale d'Emanuel. 4159 4159 IF (cvflag_ice) THEN 4160 4160 … … 4170 4170 4171 4171 !JAM Attention: evap=sigt*E 4172 ! Modification: evap devient l' évaporation en milieu de couche4173 ! car n écessaire dans cv3_yield4174 ! Du coup, il faut modifier pas mal d' équations...4172 ! Modification: evap devient l'evaporation en milieu de couche 4173 ! car necessaire dans cv3_yield 4174 ! Du coup, il faut modifier pas mal d'equations... 4175 4175 ! et l'expression de afac qui devient afac1 4176 4176 ! revap=sqrt((prec(i+1)+prec(i))/2) … … 4191 4191 !JYG Dans sa formulation originale, Emanuel calcule l'evaporation par: 4192 4192 ! c evap(il,i)=sigt*afac*revap 4193 ! ce qui n'est pas correct. Dans cv_routines, la formulation a étémodifiee.4193 ! ce qui n'est pas correct. Dans cv_routines, la formulation a ete modifiee. 4194 4194 ! Ici,l'evaporation evap est simplement calculee par l'equation de 4195 4195 ! conservation. … … 4525 4525 #ifdef ISO 4526 4526 #ifdef ISOVERIF 4527 ! verif des inputs àappel stewart4527 ! verif des inputs a appel stewart 4528 4528 do il=1,ncum 4529 4529 if (i.le.inb(il) .and. lwork(il)) then … … 4543 4543 enddo 4544 4544 #endif 4545 ! appel de appel_stewart_vectoris é4545 ! appel de appel_stewart_vectorise 4546 4546 call appel_stewart_vectall_np(lwork,ncum, & 4547 4547 & ph,t,evap,xtwdtrain, & … … 4611 4611 #endif 4612 4612 4613 ! équivalent isotopique de rp(il,i)=amin1(rp(il,i),rs(il,i))4613 ! equivalent isotopique de rp(il,i)=amin1(rp(il,i),rs(il,i)) 4614 4614 do il=1,ncum 4615 4615 if (i.lt.inb(il) .and. lwork(il)) then … … 4651 4651 #endif 4652 4652 rpprec(il,i)=rs(il,i) 4653 ! sous cas rajout éle 11dec 2011. Normalement, pas utile4653 ! sous cas rajoute le 11dec 2011. Normalement, pas utile 4654 4654 else if (rp(il,i).eq.0.0) then 4655 4655 do ixt=1,ntraciso … … 4864 4864 real xtbx(ntraciso), xtawat(ntraciso,nloc) 4865 4865 ! cam debug 4866 ! pour l'homog énéisation sous le nuage:4866 ! pour l'homogeneisation sous le nuage: 4867 4867 real bxtsum(ntraciso,nloc), fxtsum(ntraciso,nloc) 4868 4868 #ifdef DIAGISO 4869 ! diagnostiques juste: tendance des diff érents processus4869 ! diagnostiques juste: tendance des differents processus 4870 4870 real fxt_detrainement(niso,nloc,nd) 4871 4871 real fxt_fluxmasse(niso,nloc,nd) … … 4917 4917 #ifdef ISO 4918 4918 ! cam debug 4919 ! write(*,*) 'cv3_routines 3082: entr ée dans cv3_yield'4919 ! write(*,*) 'cv3_routines 3082: entree dans cv3_yield' 4920 4920 ! en cam debug 4921 4921 do ixt = 1, ntraciso … … 4994 4994 END DO 4995 4995 #ifdef ISO 4996 ! on initialise mieux fr et fxt par securit é4996 ! on initialise mieux fr et fxt par securite 4997 4997 fr(:,:)=0.0 4998 4998 fxt(:,:,:)=0.0 … … 5845 5845 else ! taggage des ddfts: 5846 5846 ! la formule pour fq_ddft suppose que le ddft est en RP. Ce n'est pas le 5847 ! cas pour le water tagging puisqu'il y a conversion des mol écules5848 ! blances entrain ées en molécule rouges.5847 ! cas pour le water tagging puisqu'il y a conversion des molecules 5848 ! blances entrainees en molecule rouges. 5849 5849 ! Il faut donc prendre en compte ce taux de conversion quand 5850 5850 ! entrainement d'env vers ddft … … 5855 5855 ! : -conversion(iiso) 5856 5856 5857 ! Pb: quand on discretise, dqp/dt n'est pas v érifée numériquement.5858 ! on se retrouve donc avec des d Ye/dt diff érents de 0 même si ye=0 ( on5859 ! note X les mol écules poubelles et Y les molécules ddfts).5857 ! Pb: quand on discretise, dqp/dt n'est pas verifee numeriquement. 5858 ! on se retrouve donc avec des d Ye/dt differents de 0 meme si ye=0 ( on 5859 ! note X les molecules poubelles et Y les molecules ddfts). 5860 5860 5861 5861 ! Solution alternative: Dans le cas entrainant, Ye ne varie que par 5862 5862 ! ascendance compensatoire des ddfts et par perte de Ye vers le ddft. On 5863 ! calcule donc ce terme directement avec sch éma amont:5864 5865 ! ajout d éjà de l'évap5863 ! calcule donc ce terme directement avec schema amont: 5864 5865 ! ajout deja de l'evap 5866 5866 do ixt = 1+niso,ntraciso 5867 5867 fxt(ixt,il,i)=fxt(ixt,il,i) & … … 5941 5941 #endif 5942 5942 else !if (abs(dXe).gt.ridicule) then 5943 ! dans ce cas, fxtXe doit être faible5943 ! dans ce cas, fxtXe doit etre faible 5944 5944 5945 5945 #ifdef ISOVERIF … … 5957 5957 fxt(ixt,il,i)=fxt(ixt,il,i)+fxtXe(iiso) 5958 5958 else !if (izone.eq.izone_poubelle) then 5959 ! pas de tendance pour ce tag l à5959 ! pas de tendance pour ce tag la 5960 5960 endif !if (izone.eq.izone_poubelle) then 5961 5961 endif !if ((izone.ne.izone_revap).and. … … 5971 5971 5972 5972 else !if (mp(il,i).gt.mp(il,i+1)) then 5973 ! cas d étrainant: pas de problèmes5973 ! cas detrainant: pas de problemes 5974 5974 do ixt=1+niso,ntraciso 5975 5975 fxt(ixt,il,i)=fxt(ixt,il,i) & … … 6176 6176 ! on change le traitement de cette ligne le 8 mai 2009: 6177 6177 ! avant, on avait: xtawat=xtelij(il,k,i)-(1.-xtep(il,i))*xtclw(il,i) 6178 ! c'est àdire que Rawat=Relij+(1-ep)*clw/awat*(Relij-Rclw)6179 ! si Relij!=Rclw, alors un fractionnement isotopique non physique était6178 ! c'est a dire que Rawat=Relij+(1-ep)*clw/awat*(Relij-Rclw) 6179 ! si Relij!=Rclw, alors un fractionnement isotopique non physique etait 6180 6180 ! introduit. 6181 ! En fait, awat repr ésente le surplus de condensat dans le mélange par6182 ! rapport àcelui restant dans la colonne adiabatique6183 ! ce surplus à la même compo que le elij, sans fractionnement.6184 ! d'o ùle nouveau traitement ci-dessous.6181 ! En fait, awat represente le surplus de condensat dans le melange par 6182 ! rapport a celui restant dans la colonne adiabatique 6183 ! ce surplus a la meme compo que le elij, sans fractionnement. 6184 ! d'ou le nouveau traitement ci-dessous. 6185 6185 if (elij(il,k,i).gt.0.0) then 6186 6186 do ixt = 1, ntraciso 6187 6187 xtawat(ixt,il)=awat(il)*(xtelij(ixt,il,k,i)/elij(il,k,i)) 6188 ! xtawat(ixt)=amax1(xtawat(ixt),0.0) ! pas n écessaire6188 ! xtawat(ixt)=amax1(xtawat(ixt),0.0) ! pas necessaire 6189 6189 enddo !do ixt = 1, ntraciso 6190 6190 else !if (elij(il,k,i).gt.0.0) then 6191 6191 ! normalement, si elij(il,k,i)<=0, alors awat=0 6192 ! on le v érifie. Si c'est vrai -> xtawat=0 aussi6192 ! on le verifie. Si c'est vrai -> xtawat=0 aussi 6193 6193 #ifdef ISOVERIF 6194 6194 call iso_verif_egalite(awat(il),0.0,'cv3_yield 3779') … … 6850 6850 fq_detrainement(il, i) = fq_detrainement(il, i)/alpha_qpos(il) 6851 6851 do ixt=1,ntraciso 6852 f q_ddft(ixt,il, i) = fq_ddft(ixt,il, i)/alpha_qpos(il)6853 f q_evapprecip(ixt,il, i) = fq_evapprecip(ixt,il, i)/alpha_qpos(il)6854 f q_fluxmasse(ixt,il, i) = fq_fluxmasse(ixt,il, i)/alpha_qpos(il)6855 f q_detrainement(ixt,il, i) = fq_detrainement(ixt,il, i)/alpha_qpos(il)6852 fxt_ddft(ixt,il, i) = fxt_ddft(ixt,il, i)/alpha_qpos(il) 6853 fxt_evapprecip(ixt,il, i) = fxt_evapprecip(ixt,il, i)/alpha_qpos(il) 6854 fxt_fluxmasse(ixt,il, i) = fxt_fluxmasse(ixt,il, i)/alpha_qpos(il) 6855 fxt_detrainement(ixt,il, i) = fxt_detrainement(ixt,il, i)/alpha_qpos(il) 6856 6856 enddo ! do ixt=1,ntraciso 6857 6857 #endif … … 7179 7179 ENDDO ! k 7180 7180 7181 ! 14/01/15 AJ delta n'a rien à faire là...7181 ! 14/01/15 AJ delta n'a rien a faire la... 7182 7182 DO il = 1, ncum ! cld 7183 7183 !! IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld … … 7195 7195 7196 7196 ! IM cf. FH 7197 ! 14/01/15 AJ ne correspond pas à ce qui a été codépar JYG et SB7197 ! 14/01/15 AJ ne correspond pas a ce qui a ete code par JYG et SB 7198 7198 7199 7199 IF (iflag_clw==0) THEN ! cld … … 7290 7290 7291 7291 ! fraction deau condensee dans les melanges convertie en precip : epm 7292 ! et eau condens ée précipitée dans masse d'air saturé: l_m*dM_m/dzdz.dzdz7292 ! et eau condensee precipitee dans masse d'air sature : l_m*dM_m/dzdz.dzdz 7293 7293 DO j = 1, nl 7294 7294 DO k = 1, nl … … 7576 7576 7577 7577 ! On fait varier epmax en fn de la cape 7578 ! Il faut donc recalculer ep, et hp qui a d éjà été calculéet7579 ! qui en d épend7580 ! Toutes les autres variables fn de ep sont calcul ées plus bas.7578 ! Il faut donc recalculer ep, et hp qui a deja ete calcule et 7579 ! qui en depend 7580 ! Toutes les autres variables fn de ep sont calculees plus bas. 7581 7581 7582 7582 include "cvthermo.h" … … 7613 7613 7614 7614 ! il faut calculer la cape: on fait un calcule simple car tant qu'on ne 7615 ! connait pas ep, on ne connait pas les m élanges, ddfts etc... qui sont7615 ! connait pas ep, on ne connait pas les melanges, ddfts etc... qui sont 7616 7616 ! necessaires au calcul de la cape dans la nouvelle physique 7617 7617 -
LMDZ6/trunk/libf/phylmdiso/infotrac_phy.F90
r4048 r4050 51 51 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: niadv ! equivalent dyn / physique 52 52 !$OMP THREADPRIVATE(niadv) 53 54 ! CRisi: tableaux de fils55 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqfils56 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqdesc ! nombres de fils + nombre de tous les petits fils sur toutes les générations57 INTEGER, SAVE :: nqdesc_tot58 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqfils59 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iqpere60 !$OMP THREADPRIVATE(nqfils,nqdesc,nqdesc_tot,iqfils,iqpere)61 53 62 54 ! conv_flg(it)=0 : convection desactivated for tracer number it … … 84 76 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqiso ! donne indice iq en fn de (ixt,phase) 85 77 !$OMP THREADPRIVATE(iqiso) 86 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_num ! donne numéro iso entre 1 et niso_possibles en fn de nqtot87 !$OMP THREADPRIVATE(iso_num)88 78 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_indnum ! donne numéro iso entre 1 et niso effectif en fn de nqtot 89 79 !$OMP THREADPRIVATE(iso_indnum) 90 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: zone_num ! donne numéro de la zone de tracage en fn de nqtot91 !$OMP THREADPRIVATE(zone_num)92 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: phase_num ! donne numéro de la zone de tracage en fn de nqtot93 !$OMP THREADPRIVATE(phase_num)94 80 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 95 81 !$OMP THREADPRIVATE(indnum_fn_num) … … 106 92 SUBROUTINE init_infotrac_phy(nqtot_,nqo_,nbtr_,nqtottr_,nqCO2_,tracers_,type_trac_,& 107 93 niadv_,conv_flg_,pbl_flg_,solsym_,& 108 nqfils_,nqdesc_,nqdesc_tot_,iqfils_,iqpere_,&109 94 ok_isotopes_,ok_iso_verif_,ok_isotrac_,& 110 95 ok_init_iso_,niso_possibles_,tnat_,& 111 alpha_ideal_,use_iso_,iqiso_,iso_num_,& 112 iso_indnum_,zone_num_,phase_num_,& 96 alpha_ideal_,use_iso_,iqiso_,iso_indnum_,& 113 97 indnum_fn_num_,index_trac_,& 114 98 niso_,ntraceurs_zone_,ntraciso_,itr_indice_& … … 143 127 CHARACTER(len=*),INTENT(IN) :: solsym_(nbtr_) 144 128 ! Isotopes: 145 INTEGER,INTENT(IN) :: nqfils_(nqtot_)146 INTEGER,INTENT(IN) :: nqdesc_(nqtot_)147 INTEGER,INTENT(IN) :: nqdesc_tot_148 INTEGER,INTENT(IN) :: iqfils_(nqtot_,nqtot_)149 INTEGER,INTENT(IN) :: iqpere_(nqtot_)150 129 LOGICAL,INTENT(IN) :: ok_isotopes_ 151 130 LOGICAL,INTENT(IN) :: ok_iso_verif_ … … 157 136 LOGICAL,INTENT(IN) :: use_iso_(niso_possibles_) 158 137 INTEGER,INTENT(IN) :: iqiso_(ntraciso_,nqo_) 159 INTEGER,INTENT(IN) :: iso_num_(nqtot_)160 138 INTEGER,INTENT(IN) :: iso_indnum_(nqtot_) 161 INTEGER,INTENT(IN) :: zone_num_(nqtot_)162 INTEGER,INTENT(IN) :: phase_num_(nqtot_)163 139 INTEGER,INTENT(IN) :: indnum_fn_num_(niso_possibles_) 164 140 INTEGER,INTENT(IN) :: index_trac_(ntraceurs_zone_,niso_) … … 216 192 217 193 IF (ok_isotopes) THEN 218 ALLOCATE(nqfils(nqtot))219 nqfils(:)=nqfils_(:)220 ALLOCATE(nqdesc(nqtot))221 nqdesc(:)=nqdesc_(:)222 nqdesc_tot=nqdesc_tot_223 ALLOCATE(iqfils(nqtot,nqtot))224 iqfils(:,:)=iqfils_(:,:)225 ALLOCATE(iqpere(nqtot))226 iqpere(:)=iqpere_(:)227 228 194 tnat(:)=tnat_(:) 229 195 alpha_ideal(:)=alpha_ideal_(:) … … 232 198 ALLOCATE(iqiso(ntraciso,nqo)) 233 199 iqiso(:,:)=iqiso_(:,:) 234 ALLOCATE(iso_num(nqtot))235 iso_num(:)=iso_num_(:)236 200 ALLOCATE(iso_indnum(nqtot)) 237 201 iso_indnum(:)=iso_indnum_(:) 238 ALLOCATE(zone_num(nqtot))239 zone_num(:)=zone_num_(:)240 ALLOCATE(phase_num(nqtot))241 phase_num(:)=phase_num_(:)242 202 243 203 indnum_fn_num(:)=indnum_fn_num_(:) … … 255 215 write(*,*) 'itr_indice=',itr_indice 256 216 #ifdef ISOVERIF 257 write(*,*) 'iso_ num=',iso_num217 write(*,*) 'iso_iName=',tracers(:)%iso_iName 258 218 #endif 259 219 -
LMDZ6/trunk/libf/phylmdiso/isotopes_verif_mod.F90
r4033 r4050 2123 2123 end function iso_verif_tracdD_choix_nostop 2124 2124 2125 subroutine iso_verif_trac17_q_deltaD(x,err_msg) 2126 use isotrac_mod, only: nzone_temp,option_traceurs 2127 USE infotrac_phy, ONLY: ntraciso 2128 implicit none 2129 2130 ! inputs 2131 real x(ntraciso) 2132 character*(*) err_msg 2133 ! local 2134 integer iso_verif_tag17_q_deltaD_chns 2135 2136 if ((option_traceurs.eq.17).or. & 2137 & (option_traceurs.eq.18)) then 2138 if (nzone_temp.ge.5) then 2139 if (iso_verif_tag17_q_deltaD_chns(x,err_msg).eq.1) then 2140 stop 2141 endif 2142 endif 2143 endif !if (option_traceurs.eq.17) then 2144 2145 end subroutine iso_verif_trac17_q_deltaD 2125 INTEGER FUNCTION iso_verif_tag17_q_deltaD_chns(x,err_msg) RESULT(res) 2126 USE infotrac_phy, ONLY: index_trac, ntraciso 2127 USE isotopes_mod, ONLY: iso_HDO, iso_eau, ridicule 2128 USE isotrac_mod, ONLY: nzone_temp, option_traceurs 2129 IMPLICIT NONE 2130 REAL, INTENT(IN) :: x(ntraciso) 2131 CHARACTER(LEN=*), INTENT(IN) :: err_msg 2132 INTEGER :: ieau, ixt, ieau1 2133 res = 0 2134 IF(ALL([17,18]/=option_traceurs)) RETURN 2135 !--- Check whether * deltaD(highest tagging layer) < 200 permil 2136 ! * q < 2137 ieau=index_trac(nzone_temp,iso_eau) 2138 ixt=index_trac(nzone_temp,iso_HDO) 2139 IF(x(ieau)>ridicule) THEN 2140 IF(iso_verif_positif_nostop(-200.0-deltaD(x(ixt)/x(ieau)), err_msg//': deltaDt05 trop fort')==1) THEN 2141 res=1; write(*,*) 'x=',x 2142 END IF 2143 END IF 2144 IF(iso_verif_positif_nostop(2.0e-3-x(ieau),err_msg//': qt05 trop fort')==1) THEN 2145 res=1; write(*,*) 'x=',x 2146 END IF 2147 !--- Check whether q is small ; then, qt01 < 10% 2148 IF(x(iso_eau)<2.0e-3) THEN 2149 ieau1= index_trac(1,iso_eau) 2150 IF(iso_verif_positif_nostop(0.1-(x(ieau1)/x(iso_eau)),err_msg//': qt01 trop abondant')==1) THEN 2151 res=1; write(*,*) 'x=',x 2152 END IF 2153 END IF 2154 END FUNCTION iso_verif_tag17_q_deltaD_chns 2155 2156 SUBROUTINE iso_verif_trac17_q_deltaD(x,err_msg) 2157 USE isotrac_mod, ONLY: nzone_temp, option_traceurs 2158 USE infotrac_phy, ONLY: ntraciso 2159 IMPLICIT NONE 2160 REAL, INTENT(IN) :: x(ntraciso) 2161 CHARACTER(LEN=*), INTENT(IN) :: err_msg 2162 IF(ALL([17,18]/=option_traceurs)) RETURN 2163 IF(nzone_temp>=5) THEN 2164 IF(iso_verif_tag17_q_deltaD_chns(x,err_msg)==1) STOP 2165 END IF 2166 END SUBROUTINE iso_verif_trac17_q_deltaD 2146 2167 2147 2168 subroutine iso_verif_traceur(x,err_msg) … … 2676 2697 2677 2698 end function iso_verif_traceur_jm_nostop 2678 2679 function iso_verif_tag17_q_deltaD_chns(x,err_msg)2680 USE infotrac_phy, ONLY: index_trac,ntraciso2681 use isotopes_mod, ONLY: iso_HDO,iso_eau,ridicule2682 use isotrac_mod, only: nzone_temp,option_traceurs2683 implicit none2684 2685 ! inputs2686 real x(ntraciso)2687 character*(*) err_msg2688 ! output2689 integer iso_verif_tag17_q_deltaD_chns2690 ! locals2691 !integer iso_verif_positif_nostop2692 !real deltaD2693 integer ieau,ixt,ieau12694 2695 iso_verif_tag17_q_deltaD_chns=02696 2697 if ((option_traceurs.eq.17).or. &2698 & (option_traceurs.eq.18)) then2699 ! verifier que deltaD du tag de la couche la plus haute <2700 ! 200 permil, et vérifier que son q est inférieur à2701 ieau=index_trac(nzone_temp,iso_eau)2702 ixt=index_trac(nzone_temp,iso_HDO)2703 2704 if (x(ieau).gt.ridicule) then2705 if (iso_verif_positif_nostop(-200.0-deltaD(x(ixt)/x(ieau)), &2706 & err_msg//': deltaDt05 trop fort').eq.1) then2707 write(*,*) 'x=',x2708 iso_verif_tag17_q_deltaD_chns=12709 endif !if (iso_verif_positif_nostop(-200.0-deltaD(x(ixt),x(ieau)),2710 endif !if (x(ieau).gt.ridicule) then2711 2712 if (iso_verif_positif_nostop(2.0e-3-x(ieau), &2713 & err_msg//': qt05 trop fort').eq.1) then2714 write(*,*) 'x=',x2715 iso_verif_tag17_q_deltaD_chns=12716 endif !if (iso_verif_positif_nostop(1.0e-3-x(ieau),2717 2718 ! on vérifie que si q est petit, alors qt01 fait moins de 10%2719 if (x(iso_eau).lt.2.0e-3) then2720 ieau1= index_trac(1,iso_eau)2721 if (iso_verif_positif_nostop(0.1-(x(ieau1)/x(iso_eau)), &2722 & err_msg//': qt01 trop abondant').eq.1) then2723 write(*,*) 'x=',x2724 iso_verif_tag17_q_deltaD_chns=12725 endif ! if (iso_verif_positif(0.1-(x(ixt)/x(ieau)),2726 endif !if (x(ieau).lt.2.0e-3) then2727 2728 endif !if (option_traceurs.eq.17) then2729 2730 end function iso_verif_tag17_q_deltaD_chns2731 2699 2732 2700 subroutine iso_verif_tag17_q_deltaD_vect(x,n,m,err_msg) -
LMDZ6/trunk/libf/phylmdiso/phyetat0.F90
r4046 r4050 45 45 use config_ocean_skin_m, only: activate_ocean_skin 46 46 #ifdef ISO 47 USE infotrac_phy, ONLY: n traciso,niso,iso_num47 USE infotrac_phy, ONLY: niso 48 48 USE isotopes_routines_mod, ONLY: phyisoetat0 49 49 USE isotopes_mod, ONLY: iso_eau -
LMDZ6/trunk/libf/phylmdiso/phys_output_mod.F90
r4048 r4050 35 35 USE iophy 36 36 USE dimphy 37 USE infotrac_phy, ONLY: nqtot, nqo, niadv, tracers, type_trac, &37 USE infotrac_phy, ONLY: nqtot, nqo, niadv, tracers, type_trac, maxlen, & 38 38 nqtottr,itr_indice ! C Risi 39 USE strings_mod, ONLY: maxlen40 39 USE ioipsl 41 40 USE phys_cal_mod, only : hour, calend … … 537 536 write(lunout,*) 'itr_indice=',itr_indice 538 537 ! IF (nqtot>=nqo+1) THEN 539 IF (nqtottr>=1) THEN540 538 ! 541 539 !DO iq=nqo+1,nqtot … … 579 577 tnam = 'cum'//TRIM(tracers(iiq)%name); o_trac_cum(itr)= ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)]) 580 578 ENDDO 581 ENDIF582 579 583 580 ENDDO ! iff -
LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
r4040 r4050 123 123 #ifdef ISO 124 124 USE infotrac_phy, ONLY: & 125 iqiso,iso_ num,iso_indnum,zone_num,ok_isotrac, &125 iqiso,iso_indnum,tracers,ok_isotrac, & 126 126 niso,ntraciso,nqtottr,itr_indice ! ajout C Risi pour isos 127 127 USE isotopes_mod, ONLY: iso_eau,iso_HDO,iso_O18,iso_O17,iso_HTO, & … … 141 141 & iso_verif_aberrant_choix,iso_verif_positif, & 142 142 & iso_verif_positif_choix_vect,iso_verif_o18_aberrant_nostop, & 143 & iso_verif_init, 143 & iso_verif_init,& 144 144 & iso_verif_positif_strict_nostop,iso_verif_O18_aberrant_enc_vect2D 145 145 #endif … … 155 155 & iso_verif_traceur_justmass,iso_verif_traceur_vect, & 156 156 & iso_verif_trac17_q_deltad,iso_verif_trac_masse_vect, & 157 & iso_verif_t racpos_choix_nostop157 & iso_verif_tag17_q_deltaD_vect, iso_verif_tracpos_choix_nostop 158 158 #endif 159 159 #endif … … 2366 2366 #endif 2367 2367 if (ixt.gt.niso) then 2368 write(*,*) 'izone,iiso=', zone_num(iqiso(ixt,ivap)),iso_indnum(iqiso(ixt,ivap))2368 write(*,*) 'izone,iiso=',tracers(iqiso(ixt,ivap))%iso_iZone,iso_indnum(iqiso(ixt,ivap)) 2369 2369 endif 2370 2370 DO k = 1, klev
Note: See TracChangeset
for help on using the changeset viewer.