Changeset 3195
- Timestamp:
- Jan 31, 2024, 4:36:51 PM (11 months ago)
- Location:
- trunk/LMDZ.PLUTO
- Files:
-
- 21 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.PLUTO/changelog.txt
r3184 r3195 1838 1838 == 25/01/2024 == AF 1839 1839 Added LMDZ.PLUTO, a copy of the generic model cleaned of some 1840 unnecessary 1840 unnecessary modules -
trunk/LMDZ.PLUTO/deftank/callphys.def
r3193 r3195 6 6 diurnal = .true. 7 7 # Seasonal cycle ? if season=false, Ls stays constant, to value set in "start" 8 season = .true. 8 season = .true. 9 9 # Tidally resonant orbit ? must have diurnal=false, correct rotation rate in newstart 10 10 tlocked = .false. … … 82 82 Fat1AU = 1366.0 83 83 84 ## Tracer and aerosol options 84 ## Tracer and aerosol options 85 85 ## ~~~~~~~~~~~~~~~~~~~~~~~~~~ 86 86 # Gravitational sedimentation of tracers (just H2O ice for now) ? 87 87 sedimentation = .true. 88 88 89 naerkind = 189 naerkind = 0 90 90 hazeprop = optprop_tholins_fractal100nm.dat 91 91 … … 105 105 ## extra non-standard definitions for Earth 106 106 ######################################################################### 107 108 ## Tracer and aerosol options 107 108 ## Tracer and aerosol options 109 109 ## ~~~~~~~~~~~~~~~~~~~~~~~~~~ 110 110 # atm mass update due to tracer evaporation/condensation? … … 123 123 124 124 125 ## Condensation 125 ## Condensation 126 126 ## ~~~~~~~~~~~ 127 127 # call N2 condensation ? -
trunk/LMDZ.PLUTO/deftank/traceur.def
r3193 r3195 1 1 2 n2 1 0 -
trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/lect_start_archive.F
r3184 r3195 7 7 8 8 USE comsoil_h, ONLY: nsoilmx, layer, mlayer, volcapa, inertiedat 9 USE tracer_h, ONLY: igcm_ n2_ice9 USE tracer_h, ONLY: igcm_haze 10 10 USE infotrac, ONLY: tname, nqtot 11 11 ! USE slab_ice_h, ONLY: noceanmx … … 39 39 40 40 41 c Variables pour les lectures des fichiers "ini" 41 c Variables pour les lectures des fichiers "ini" 42 42 c-------------------------------------------------- 43 43 ! INTEGER sizei, … … 53 53 ! character (len=50) :: tmpname 54 54 55 c Variable histoire 55 c Variable histoire 56 56 c------------------ 57 57 REAL,INTENT(OUT) :: vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants … … 59 59 REAL,INTENT(OUT) :: q(iip1,jjp1,llm,nqtot) 60 60 61 c Physique sur grille scalaire 61 c Physique sur grille scalaire 62 62 c---------------------------- 63 63 … … 112 112 113 113 114 c Variable de l'ancienne grille 114 c Variable de l'ancienne grille 115 115 c--------------------------------------------------------- 116 116 … … 150 150 logical :: therminertia_3D=.true. ! flag 151 151 ! therminertia_3D=.true. if thermal inertia is 3D and read from datafile 152 c Variable intermediaires iutilise pour l'extrapolation verticale 152 c Variable intermediaires iutilise pour l'extrapolation verticale 153 153 c---------------------------------------------------------------- 154 real, dimension(:,:,:), allocatable :: var,varp1 154 real, dimension(:,:,:), allocatable :: var,varp1 155 155 real, dimension(:), allocatable :: oldgrid, oldval 156 156 real, dimension(:), allocatable :: newval … … 158 158 real,intent(out) :: surfith(iip1,jjp1) ! surface thermal inertia 159 159 ! surface thermal inertia at old horizontal grid resolution 160 real, dimension(:,:), allocatable :: surfithold 160 real, dimension(:,:), allocatable :: surfithold 161 161 162 162 character(len=30) :: txt ! to store some text … … 172 172 !----------------------------------------------------------------------- 173 173 174 ! 1.2 Read the various dimension lengths of data in file 174 ! 1.2 Read the various dimension lengths of data in file 175 175 176 176 ierr= NF_INQ_DIMID(nid,"Time",dimid) … … 245 245 246 246 ! 1.3 Report dimensions 247 247 248 248 write(*,*) "Start_archive dimensions:" 249 249 write(*,*) "longitude: ",imold … … 259 259 endif 260 260 write(*,*) "time lenght: ",timelen 261 write(*,*) 261 write(*,*) 262 262 263 263 !----------------------------------------------------------------------- … … 311 311 312 312 C----------------------------------------------------------------------- 313 c 3.1. Lecture du tableau des parametres du run 313 c 3.1. Lecture du tableau des parametres du run 314 314 c (pour la lecture ulterieure de "ptotalold" et "n2icetotalold") 315 315 c----------------------------------------------------------------------- … … 355 355 PRINT*, "lect_start_archive: Field <rlatu> not found" 356 356 CALL abort 357 ENDIF 357 ENDIF 358 358 #ifdef NC_DOUBLE 359 359 ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold) … … 440 440 c 3.4 Read Soil layers depths 441 441 c----------------------------------------------------------------------- 442 442 443 443 ierr=NF_INQ_VARID(nid,"soildepth",nvarid) 444 444 if (ierr.ne.NF_NOERR) then … … 540 540 ptotalold = tab_cntrl(tab0+49) 541 541 n2icetotalold = tab_cntrl(tab0+50) 542 542 543 543 c----------------------------------------------------------------------- 544 544 c 4. Lecture du temps et choix 545 545 c----------------------------------------------------------------------- 546 546 547 547 c lecture du temps 548 548 c … … 559 559 IF (ierr .NE. NF_NOERR) THEN 560 560 ierr = NF_INQ_VARID (nid, "temps", nvarid) 561 endif 561 endif 562 562 #ifdef NC_DOUBLE 563 563 ierr = NF_GET_VAR_DOUBLE(nid, nvarid, timelist) … … 584 584 585 585 6 FORMAT(i7,i7,f9.3) 586 586 587 587 write(*,*) 588 588 write(*,*) 'Choice for the date' … … 595 595 endif 596 596 end do 597 597 598 598 if (memo.eq.0) then 599 599 write(*,*) … … 618 618 c 5.1 Lecture des champs 2D (n2ice, emis,ps,tsurf,Tg[10], qsurf) 619 619 c----------------------------------------------------------------------- 620 620 621 621 start=(/1,1,memo,0/) 622 622 count=(/imold+1,jmold+1,1,0/) 623 623 624 624 ierr = NF_INQ_VARID (nid, "emis", nvarid) 625 625 IF (ierr .NE. NF_NOERR) THEN … … 756 756 ! PRINT*, "lect_start_archive: Failed loading <sea_ice>" 757 757 ! ENDIF 758 758 759 759 ! ENDIF! ok_slab_ocean 760 760 c … … 763 763 write(*,*) 764 764 765 ! Surface tracers: 765 ! Surface tracers: 766 766 ! initialize all surface tracers to zero 767 767 qsurfold(1:imold+1,1:jmold+1,1:nqtot)=0 … … 777 777 endif 778 778 779 779 780 780 write(*,*) "lect_start_archive: loading tracer ",trim(txt) 781 781 ierr = NF_INQ_VARID (nid,txt,nvarid) … … 838 838 c 839 839 enddo ! of do isoil=1,nsoilold 840 840 841 841 ! reset 'start' and 'count' to "3D" behaviour 842 842 start=(/1,1,1,memo/) 843 843 count=(/imold+1,jmold+1,nsoilold,1/) 844 844 845 845 else 846 846 write(*,*) "lect_start_archive: loading tsoil " … … 856 856 #endif 857 857 endif ! of if (ierr.ne.NF_NOERR) 858 858 859 859 endif ! of if (olddepthdef) 860 860 … … 928 928 c 929 929 930 ! Tracers: 930 ! Tracers: 931 931 qold(1:imold+1,1:jmold+1,1:lmold,1:nqtot)=0. 932 932 … … 1040 1040 c INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables 1041 1041 c======================================================================= 1042 c Interpolation horizontale puis passage dans la grille physique pour 1043 c les variables physique 1042 c Interpolation horizontale puis passage dans la grille physique pour 1043 c les variables physique 1044 1044 c Interpolation verticale puis horizontale pour chaque variable 3D 1045 1045 c======================================================================= … … 1048 1048 c 6.1 Variable 2d : 1049 1049 c----------------------------------------------------------------------- 1050 c Relief 1050 c Relief 1051 1051 call interp_horiz (phisold,phisold_newgrid,imold,jmold,iim,jjm,1, 1052 1052 & rlonuold,rlatvold,rlonu,rlatv) 1053 1053 1054 ! N2 ice is now in qsurf(igcm_ n2_ice)1054 ! N2 ice is now in qsurf(igcm_haze) 1055 1055 ! call interp_horiz (n2iceold,n2ices,imold,jmold,iim,jjm,1, 1056 1056 ! & rlonuold,rlatvold,rlonu,rlatv) … … 1088 1088 c On assure la conservation de la masse de l'atmosphere + calottes 1089 1089 c----------------------------------------------------------------------- 1090 !AF: TODO: mass conservation: check this. haze? 1090 1091 1091 1092 ptotal = 0. … … 1096 1097 END DO 1097 1098 n2icetotal = 0. 1098 if (igcm_ n2_ice.ne.0) then1099 if (igcm_haze.ne.0) then 1099 1100 ! recast surface N2 ice on new grid 1100 call interp_horiz(qsurfold(1,1,igcm_ n2_ice),1101 & qsurfs(1,1,igcm_ n2_ice),1101 call interp_horiz(qsurfold(1,1,igcm_haze), 1102 & qsurfs(1,1,igcm_haze), 1102 1103 & imold,jmold,iim,jjm,1, 1103 1104 & rlonuold,rlatvold,rlonu,rlatv) … … 1105 1106 DO i=1,iim 1106 1107 !n2icetotal = n2icetotal + n2iceS(i,j)*aire(i,j) 1107 n2icetotal=n2icetotal+qsurfS(i,j,igcm_ n2_ice)*aire(i,j)1108 n2icetotal=n2icetotal+qsurfS(i,j,igcm_haze)*aire(i,j) 1108 1109 END DO 1109 1110 END DO … … 1115 1116 write(*,*)'Old grid: atmospheric mass :',ptotalold 1116 1117 write(*,*)'New grid: atmospheric mass :',ptotal 1117 write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold 1118 write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold 1118 1119 write(*,*) 1119 1120 write(*,*)'Old grid: mass of N2 ice:',n2icetotalold … … 1131 1132 END DO 1132 1133 1133 if ( n2icetotalold.gt.0.) then 1134 if ( n2icetotalold.gt.0.) then 1134 1135 ! DO j=1,jjp1 1135 1136 ! DO i=1,iip1 … … 1157 1158 &f soil thermal inertia; might be wiser to reset it.' 1158 1159 write(*,*) 1159 1160 1160 1161 do i=1,imold+1 1161 1162 do j=1,jmold+1 … … 1177 1178 ! We have inertiedatold 1178 1179 if((imold.ne.iim).or.(jmold.ne.jjm)) then 1179 write(*,*)'lect_start_archive: WARNING: horizontal interpolation 1180 write(*,*)'lect_start_archive: WARNING: horizontal interpolation 1180 1181 &of thermal inertia; might be better to reset it.' 1181 1182 write(*,*) 1182 1183 endif 1183 1184 1184 1185 ! Do horizontal interpolation 1185 1186 if (depthinterpol) then … … 1210 1211 call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid, 1211 1212 & inertiedatS,inertiedat) 1212 1213 1213 1214 c----------------------------------------------------------------------- 1214 1215 c 6.2.2 Soil temperature … … 1282 1283 deallocate(oldval) 1283 1284 deallocate(newval) 1284 1285 1285 1286 else 1286 1287 tsoiloldnew(:,:,:)=tsoilold(:,:,:) … … 1303 1304 c 6.4 Variable 3d : 1304 1305 c----------------------------------------------------------------------- 1305 1306 1306 1307 c temperatures atmospheriques 1307 1308 write (*,*) 'lect_start_archive: told ', told (1,jmold+1,1) ! INFO … … 1321 1322 & rlonuold,rlatvold,rlonu,rlatv) 1322 1323 call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,du_nonoro_gwdS,du_nonoro_gwd) 1323 1324 1324 1325 call interp_vert 1325 1326 & (dv_nonoro_gwdold,var,lmold,llm,apsold,bpsold,aps,bps, … … 1328 1329 & rlonuold,rlatvold,rlonu,rlatv) 1329 1330 call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,dv_nonoro_gwdS,dv_nonoro_gwd) 1330 1331 1331 1332 call interp_vert 1332 1333 & (east_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps, … … 1335 1336 & rlonuold,rlatvold,rlonu,rlatv) 1336 1337 call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,east_gwstressS,east_gwstress) 1337 1338 1338 1339 call interp_vert 1339 1340 & (west_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps, … … 1379 1380 end do 1380 1381 end do 1381 end do 1382 end do 1382 1383 write (*,*) 'lect_start_archive: ucov ', ucov (1,2,1) ! INFO 1383 1384 c write(48,*) 'ucov',ucov … … 1410 1411 & rlonuold,rlatvold,rlonu,rlatv) 1411 1412 enddo 1412 cccccccccccccccccccccccccccccc 1413 c make sure that sum of q = 1 1414 c dominent species is = 1 - sum(all other species) 1415 cccccccccccccccccccccccccccccc 1413 cccccccccccccccccccccccccccccc 1414 c make sure that sum of q = 1 1415 c dominent species is = 1 - sum(all other species) 1416 cccccccccccccccccccccccccccccc 1416 1417 c iqmax=1 1417 c 1418 c 1418 1419 c if (nqold.gt.10) then 1419 1420 c do l=1,llm … … 1428 1429 c qtot(i,j,l)=0 1429 1430 c do iq=1,nqold 1430 c if (iq.ne.iqmax) then 1431 c q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq) 1431 c if (iq.ne.iqmax) then 1432 c q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq) 1432 1433 c endif 1433 1434 c enddo !iq … … 1437 1438 c $ qtot(i,j,l) 1438 1439 c enddo !iq 1439 c enddo !i 1440 c enddo !j 1441 c enddo !l 1440 c enddo !i 1441 c enddo !j 1442 c enddo !l 1442 1443 c endif 1443 1444 ccccccccccccccccccccccccccccccc … … 1451 1452 end do 1452 1453 enddo 1453 1454 1454 1455 ! call gr_dyn_fi (1,iim+1,jjm+1,ngrid,n2ices,n2ice) 1455 ! no need to transfer "n2ice" any more; it is in qsurf(igcm_ n2_ice)1456 ! no need to transfer "n2ice" any more; it is in qsurf(igcm_haze) 1456 1457 1457 1458 endif !! if nqtot .ne. 0 -
trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90
r3184 r3195 10 10 11 11 use radinc_h, only : L_TAUMAX,naerkind 12 use aerosol_mod, only: iaero_nlay, iaero_generic, & 13 iaero_aurora, iaero_back2lay, iaero_n2, & 14 iaero_dust, iaero_h2so4, & 15 iaero_nh3, i_rgcs_ice, noaero 12 use aerosol_mod, only: iaero_haze, i_haze, iaero_generic, & 13 noaero 16 14 USE tracer_h, only: noms,rho_n2,rho_ice,rho_q,mmol 17 15 use comcstfi_mod, only: g, pi, mugaz, avocado 18 16 use geometry_mod, only: latitude 19 17 use callkeys_mod, only: aerofixn2,kastprof,cloudlvl, & 20 dusttau, &21 18 pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo, & 22 19 pres_bottom_strato,pres_top_strato,obs_tau_col_strato, & 23 tau_nh3_cloud, pres_nh3_cloud, & 24 nlayaero, aeronlay_tauref, aeronlay_choice, & 20 nlayaero, aeronlay_tauref, aeronlay_choice, & 25 21 aeronlay_pbot, aeronlay_ptop, aeronlay_sclhght, & 26 22 aerogeneric … … 29 25 30 26 !================================================================== 31 ! 27 ! 32 28 ! Purpose 33 29 ! ------- 34 30 ! Compute aerosol optical depth in each gridbox. 35 ! 31 ! 36 32 ! Authors 37 ! ------- 33 ! ------- 38 34 ! F. Forget 39 ! F. Montmessin (water ice scheme) 35 ! F. Montmessin (water ice scheme) 40 36 ! update J.-B. Madeleine (2008) 41 37 ! dust removal, simplification by Robin Wordsworth (2009) … … 44 40 ! 45 41 ! Input 46 ! ----- 42 ! ----- 47 43 ! ngrid Number of horizontal gridpoints 48 44 ! nlayer Number of layers … … 97 93 CHARACTER(LEN=20) :: tracername ! to temporarily store text 98 94 99 ! for fixed dust profiles100 real topdust, expfactor, zp101 REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling102 REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling103 104 95 real CLFtot 105 96 integer igen_ice,igen_vap ! to store the index of generic tracer … … 131 122 endif 132 123 133 if ((iaero_n2.ne.0).and.(.not.noaero)) then134 print*, 'iaero_n2= ',iaero_n2135 endif136 if (iaero_dust.ne.0) then137 print*,'iaero_dust= ',iaero_dust138 endif139 if (iaero_h2so4.ne.0) then140 print*,'iaero_h2so4= ',iaero_h2so4141 endif142 if (iaero_back2lay.ne.0) then143 print*,'iaero_back2lay= ',iaero_back2lay144 endif145 if (iaero_nh3.ne.0) then146 print*,'iaero_nh3= ',iaero_nh3147 endif148 if (iaero_nlay(1).ne.0) then149 print*,'iaero_nlay= ',iaero_nlay(:)150 endif151 if (iaero_aurora.ne.0) then152 print*,'iaero_aurora= ',iaero_aurora153 endif154 if (aerogeneric .ne. 0) then155 print*,"iaero_generic= ",iaero_generic(:)156 endif157 124 firstcall=.false. 158 125 ENDIF ! of IF (firstcall) … … 161 128 ! --------------------------------------------------------- 162 129 !================================================================== 163 ! N2 ice aerosols130 ! Haze aerosols 164 131 !================================================================== 165 132 166 if (iaero_ n2.ne.0) then167 iaer=iaero_n2133 if (iaero_haze.ne.0) then 134 iaer=iaero_haze 168 135 ! 1. Initialization 169 136 aerosol(1:ngrid,1:nlayer,iaer)=0.0 170 137 ! 2. Opacity calculation 171 if (noaero) then ! aerosol set to zero 172 aerosol(1:ngrid,1:nlayer,iaer)=0.0 173 elseif (aerofixn2.or.(i_n2ice.eq.0)) then ! N2 ice cloud prescribed 174 aerosol(1:ngrid,1:nlayer,iaer)=1.e-9 175 !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option 176 else 177 DO ig=1, ngrid 178 DO l=1,nlayer-1 ! to stop the rad tran bug 179 180 aerosol0 = & 181 ( 0.75 * QREFvis3d(ig,l,iaer) / & 182 ( rho_n2 * reffrad(ig,l,iaer) ) ) * & 183 ( pq(ig,l,i_n2ice) + 1.E-9 ) * & 184 ( pplev(ig,l) - pplev(ig,l+1) ) / g 185 aerosol0 = max(aerosol0,1.e-9) 186 aerosol0 = min(aerosol0,L_TAUMAX) 187 aerosol(ig,l,iaer) = aerosol0 188 ! aerosol(ig,l,iaer) = 0.0 189 ! print*, aerosol(ig,l,iaer) 190 ! using cloud fraction 191 ! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF)) 192 ! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX) 193 194 195 ENDDO 138 DO ig=1, ngrid 139 DO l=1,nlayer-1 ! to stop the rad tran bug 140 ! if fractal, radius doit etre equivalent sphere radius 141 aerosol0 = & 142 ( 0.75 * QREFvis3d(ig,l,iaer) / & 143 ( rho_q(i_haze) * reffrad(ig,l,iaer) ) ) * & 144 ( pq(ig,l,i_haze) + 1.E-10 ) * & 145 ( pplev(ig,l) - pplev(ig,l+1) ) / g 146 aerosol0 = max(aerosol0,1.e-10) 147 aerosol0 = min(aerosol0,L_TAUMAX) 148 aerosol(ig,l,iaer) = aerosol0 196 149 ENDDO 197 end if ! if fixed or varying198 end if ! if N2 aerosols199 !==================================================================200 ! Water ice / liquid (AF24: removed)201 !==================================================================202 203 !==================================================================204 ! Dust (AF24: removed)205 !==================================================================206 207 !==================================================================208 ! H2SO4 (AF24: removed)209 !==================================================================210 !==================================================================211 ! Two-layer aerosols (unknown composition)212 ! S. Guerlet (2013) - Modif by J. Vatant d'Ollone (2020)213 !214 ! This scheme is deprecated and left for retrocompatibility215 ! You should use the n-layer scheme below !216 !217 !==================================================================218 219 !==================================================================220 ! Saturn/Jupiter ammonia cloud = thin cloud (scale height 0.2 hard coded...)221 ! S. Guerlet (2013)222 ! JVO 20 : You should now use the generic n-layer scheme below223 !==================================================================224 225 !=========================================================================================================226 ! Generic N-layers aerosols/clouds227 ! Author : J. Vatant d'Ollone (2020)228 !229 ! Purpose: Replaces and extents the former buggy 2-layer scheme as well as hard-coded NH3 cloud230 !231 ! + Each layer can have different optical properties, size of particle ...232 ! + Enables up to n=4 layers as we apparently cannot run with more scatterers (could be worth checking...)233 ! + You have different choices for vertical profile of the aerosol layers :234 ! * aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.235 ! * aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).236 ! In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.237 ! * aeronlay_choice = ... feel free to add more cases !238 ! + Layers can overlap if needed (if you want a 'transition layer' as in the 2-scheme, just add it)239 !240 !=========================================================================================================241 242 if (iaero_nlay(1) .ne.0) then243 244 DO ia=1,nlayaero245 iaer=iaero_nlay(ia)246 247 ! a. Initialization248 aerosol(1:ngrid,1:nlayer,iaer)=0.D0249 250 ! b. Opacity calculation251 252 ! Case 1 : Follows atmospheric scale height between boundaries pressures253 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~254 IF (aeronlay_choice(ia).eq.1) THEN255 256 dp_layer(:)=0.D0257 DO ig=1,ngrid258 DO l=1,nlayer-1259 !! i. Opacity follows scale height260 IF ( pplev(ig,l).le.aeronlay_pbot(ia) .AND. &261 pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN262 aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )263 dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)264 !! ii. Outside aerosol layer boundaries: no aerosols265 ELSE266 aerosol(ig,l,iaer) = 0.D0267 ENDIF268 ENDDO269 ENDDO270 ! iii. Re-normalize to required total opacity271 DO ig=1,ngrid272 DO l=1,nlayer-1273 IF ( pplev(ig,l).le.aeronlay_pbot(ia) .AND. &274 pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN275 aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &276 * aeronlay_tauref(ia)277 ENDIF278 ENDDO279 ENDDO280 281 ! Case 2 : Follows input scale height282 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~283 ELSE IF (aeronlay_choice(ia).eq.2) THEN284 285 cloud_slope = 1.D0/aeronlay_sclhght(ia)286 pcloud_deck = aeronlay_pbot(ia)287 dp_layer(:) = 0.D0288 289 DO ig=1,ngrid290 DO l=1,nlayer-1291 !! i. Below cloud layer: no opacity292 IF (pplev(ig,l) .gt. pcloud_deck) THEN293 aerosol(ig,l,iaer) = 0.D0294 !! ii. Follows scale height above cloud deck295 ELSEIF (pplev(ig,l) .le. pcloud_deck) THEN296 aerosol(ig,l,iaer) = ((pplev(ig,l)/pcloud_deck)**(cloud_slope))297 dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)298 ENDIF299 ENDDO300 ENDDO301 ! iii. Re-normalize to required total opacity302 DO ig=1,ngrid303 DO l=1,nlayer-1304 IF (pplev(ig,l) .le. pcloud_deck) THEN305 aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &306 * aeronlay_tauref(ia)307 ENDIF308 ENDDO309 ENDDO310 311 ENDIF ! aeronlay_choice312 313 ENDDO ! loop on n aerosol layers314 315 end if ! if N-layer aerosols316 317 !==================================================================318 ! Jovian auroral aerosols (unknown composition) NON-GENERIC: vertical and meridional profile tuned to observations319 ! S. Guerlet (2015)320 !==================================================================321 322 323 if (iaero_aurora .ne.0) then324 iaer=iaero_aurora325 ! 1. Initialization326 aerosol(1:ngrid,1:nlayer,iaer)=0.D0327 pm = 2000. !!case study: maxi aerosols at 20 hPa328 ! 2. Opacity calculation329 DO ig=1,ngrid330 331 !! Test Jupiter (based on Zhang et al 2013 observations, but a bit different), decembre 2015332 DO l=1,nlayer333 aerosol(ig,l,iaer) = (pplev(ig,l)/pm)**2 * exp(-(pplev(ig,l)/pm)**2)334 ENDDO335 ENDDO336 337 ! 3. Meridional distribution, and re-normalize to observed total column338 dp_layer(:)=0.D0339 DO ig=1,ngrid340 !!Jupiter341 !!Hem sud:342 IF (latitude(ig)*180.D0/pi .lt. -45.D0 .and. latitude(ig)*180.D0/pi .gt. -70.) THEN343 obs_tau_col_aurora= 10.D0**(-0.06D0*latitude(ig)*180.D0/pi-3.4D0)344 ELSEIF (latitude(ig)*180.D0/pi .lt. -37.D0 .and. latitude(ig)*180.D0/pi .ge. -45.) THEN345 obs_tau_col_aurora= 10.D0**(-0.3D0*latitude(ig)*180.D0/pi-14.3D0)346 ELSEIF (latitude(ig)*180./pi .le. -70. ) THEN347 obs_tau_col_aurora= 10**(0.06*70.-3.4D0)348 !!Hem Nord:349 ELSEIF (latitude(ig)*180.D0/pi .gt. 30.D0 .and. latitude(ig)*180.D0/pi .lt. 70.) THEN350 obs_tau_col_aurora= 10.D0**(0.03D0*latitude(ig)*180.D0/pi-1.17D0)351 ELSEIF (latitude(ig)*180.D0/pi .gt. 22.D0 .and. latitude(ig)*180.D0/pi .le. 30.) THEN352 obs_tau_col_aurora= 10.D0**(0.3D0*latitude(ig)*180.D0/pi-9.4D0)353 ELSEIF (latitude(ig)*180.D0/pi .ge. 70.) THEN354 obs_tau_col_aurora= 10**(0.03*70.-1.17D0)355 ELSEIF (latitude(ig)*180.D0/pi .ge. -37. .and. latitude(ig)*180.D0/pi .le. 22.) THEN356 obs_tau_col_aurora = 0.001D0 !!Jupiter: mini pas a zero357 ENDIF358 359 DO l=1,nlayer360 dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora361 ENDDO362 150 ENDDO 363 364 DO ig=1,ngrid 365 DO l=1,nlayer-1 366 aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/dp_layer(ig) 367 ENDDO 368 ENDDO 369 370 371 end if ! if Auroral aerosols 372 !=========================================================================== 373 ! Radiative Generic Condensable aerosols scheme 374 ! Only used when we give aerogeneric != 0 in callphys.def 375 ! Computes the generic aerosols' opacity in the same fashion as water of 376 ! dust, using the QREFvis3d of the concerned specie 377 ! Lucas Teinturier (2022) 378 !=========================================================================== 379 if (aerogeneric .ne. 0) then ! we enter the scheme 380 do ia=1,aerogeneric 381 iaer = iaero_generic(ia) 382 ! Initialization 383 aerosol(1:ngrid,1:nlayer,iaer) = 0.D0 384 igen_ice = i_rgcs_ice(ia) 385 ! Let's loop on the horizontal and vertical grid 386 do ig=1,ngrid 387 do l=1,nlayer 388 aerosol(ig,l,iaer) = ( 0.75*QREFvis3d(ig,l,iaer) / & 389 (rho_q(igen_ice) * reffrad(ig,l,iaer)) ) * & 390 (pq(ig,l,igen_ice)+1E-9 ) * & 391 (pplev(ig,l) - pplev(ig,l+1)) /g 392 enddo !l=1,nlayer 393 enddo !ig=1,ngrid 394 enddo !ia=1,aerogeneric 395 endif !aerogeneric .ne. 0 151 !QREF est le meme dans toute la colonne par def si size uniforme 152 !print*, 'TB17: QREFvis3d=',QREFvis3d(1,:,1) 153 !print*, 'TB17: rho_q=',rho_q(i_haze) 154 !print*, 'TB17: reffrad=',reffrad(1,:,1) 155 !print*, 'TB17: pq=',pq(1,:,i_haze) 156 !print*, 'TB17: deltap=',pplev(1,1) - pplev(1,nlayer) 157 end if ! if haze aerosols 396 158 397 159 ! -------------------------------------------------------------------------- 398 160 ! Column integrated visible optical depth in each point (used for diagnostic) 399 161 400 tau_col(:)=0.0 401 do iaer = 1, naerkind 402 do l=1,nlayer 403 do ig=1,ngrid 404 tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer) 405 end do 162 tau_col(:)=0.0 163 do iaer = 1, naerkind 164 do l=1,nlayer 165 do ig=1,ngrid 166 tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer) 406 167 end do 407 168 end do 169 end do 408 170 409 !do ig=1,ngrid410 !do l=1,nlayer411 !do iaer = 1, naerkind412 !if(aerosol(ig,l,iaer).gt.1.e3)then413 !print*,'WARNING: aerosol=',aerosol(ig,l,iaer)414 !print*,'at ig=',ig,', l=',l,', iaer=',iaer415 !print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)416 !print*,'reffrad=',reffrad(ig,l,iaer)417 !endif418 !end do419 !end do420 !end do171 do ig=1,ngrid 172 do l=1,nlayer 173 do iaer = 1, naerkind 174 if(aerosol(ig,l,iaer).gt.1.e3)then 175 print*,'WARNING: aerosol=',aerosol(ig,l,iaer) 176 print*,'at ig=',ig,', l=',l,', iaer=',iaer 177 print*,'QREFvis3d=',QREFvis3d(ig,l,iaer) 178 print*,'reffrad=',reffrad(ig,l,iaer) 179 endif 180 end do 181 end do 182 end do 421 183 422 ! do ig=1,ngrid 423 ! if(tau_col(ig).gt.1.e3)then 424 ! print*,'WARNING: tau_col=',tau_col(ig) 425 ! print*,'at ig=',ig 426 ! print*,'aerosol=',aerosol(ig,:,:) 427 ! print*,'QREFvis3d=',QREFvis3d(ig,:,:) 428 ! print*,'reffrad=',reffrad(ig,:,:) 429 ! endif 430 ! end do 184 do ig=1,ngrid 185 if(tau_col(ig).gt.1.e3)then 186 print*,'WARNING: tau_col=',tau_col(ig) 187 print*,'at ig=',ig 188 print*,'aerosol=',aerosol(ig,:,:) 189 print*,'QREFvis3d=',QREFvis3d(ig,:,:) 190 print*,'reffrad=',reffrad(ig,:,:) 191 endif 192 end do 193 end subroutine aeropacity 431 194 432 end subroutine aeropacity433 434 195 end module aeropacity_mod -
trunk/LMDZ.PLUTO/libf/phypluto/aerosol_mod.F90
r3193 r3195 8 8 ! corresponding aerosol was not activated in callphys.def 9 9 ! -- otherwise a value is set via iniaerosol 10 integer, save, protected :: iaero_n2 = 0 11 integer, save, protected :: iaero_dust = 0 12 integer, save, protected :: iaero_h2so4 = 0 10 integer, save :: iaero_haze = 0 11 integer, save :: i_haze = 0 13 12 logical, save, protected :: noaero = .false. 14 !$OMP THREADPRIVATE(iaero_ n2,iaero_h2o,iaero_dust,iaero_h2so4,noaero)13 !$OMP THREADPRIVATE(iaero_haze,i_haze,noaero) 15 14 16 15 ! two-layer simple aerosol model 17 16 integer, save, protected :: iaero_back2lay = 0 18 ! NH3 cloud19 integer, save, protected :: iaero_nh3 = 020 17 ! N-layer aerosol model (replaces the 2-layer and hard-coded clouds) 21 18 integer,dimension(:), allocatable, save, protected :: iaero_nlay 22 ! Auroral aerosols 23 integer, save, protected :: iaero_aurora = 0 24 !$OMP THREADPRIVATE(iaero_back2lay,iaero_nh3,iaero_nlay,iaero_aurora) 19 !$OMP THREADPRIVATE(iaero_back2lay,iaero_nlay) 25 20 26 21 ! Generic aerosols … … 33 28 contains 34 29 35 SUBROUTINE iniaerosol 30 !================================================================== 31 subroutine haze_prof(ngrid,nlayer,zzlay,pplay,pt,reffrad,profmmr) 32 !================================================================== 33 ! Purpose 34 ! ------- 35 ! Get fixed haze properties 36 ! profile of haze (from txt file) and fixed radius profile 37 ! 38 !================================================================== 39 use radinc_h, only: naerkind 40 use datafile_mod 41 use tracer_h 42 use comcstfi_mod, only: r, pi 36 43 37 use mod_phys_lmdz_para, only : is_master 38 use radinc_h, only: naerkind 39 use tracer_h, only: n_rgcs, nqtot, is_rgcs 40 use callkeys_mod, only: aeron2, dusttau, & 41 aeroback2lay, aeronh3, nlayaero, aeronlay, & 42 aeroaurora, aerogeneric 44 !----------------------------------------------------------------------- 45 ! Arguments 46 Implicit none 43 47 44 IMPLICIT NONE 45 !======================================================================= 46 ! subject: 47 ! -------- 48 ! Initialization related to aerosols 49 ! (N2 aerosols, dust, water, chemical species, ice...) 50 ! 51 ! author: Laura Kerber, S. Guerlet 52 ! ------ 53 ! 54 !======================================================================= 48 integer,intent(in) :: ngrid 49 integer,intent(in) :: nlayer 50 real,intent(in) :: zzlay(ngrid,nlayer) 51 real,intent(in) :: pplay(ngrid,nlayer) 52 real,intent(in) :: pt(ngrid,nlayer) 53 real, intent(in) :: reffrad(ngrid,nlayer,naerkind) ! haze particles radii (m) 55 54 56 integer :: i, ia, iq55 real, intent(out) :: profmmr(ngrid,nlayer) ! mmr haze kg/kg 57 56 58 ! Special case, dyn. allocation for n-layer depending on callphys.def 59 IF(.NOT.ALLOCATED(iaero_nlay)) ALLOCATE(iaero_nlay(nlayaero)) 60 iaero_nlay(:) = 0 61 ! Do the same for iaero_generic and i_rgcs_ice 62 IF (.not. allocated(iaero_generic)) allocate(iaero_generic(aerogeneric)) 63 if (.not. allocated(i_rgcs_ice)) allocate(i_rgcs_ice(aerogeneric)) 57 ! Local variables 58 integer :: iaer,l,ig,ifine 64 59 65 ! Init of i_rgcs_ice 66 i_rgcs_ice(:) =0 67 ia = 1 68 do iq=1,nqtot 69 if (is_rgcs(iq) .eq. 1) then 70 i_rgcs_ice(ia)=iq 71 ia = ia+1 72 endif 73 enddo 60 LOGICAL firstcall 61 SAVE firstcall 62 DATA firstcall/.true./ 74 63 75 iaero_generic(:)=0 76 ia=0 77 if (aeron2) then 78 ia=ia+1 79 iaero_n2=ia 80 endif 81 if (is_master) write(*,*) '--- N2 aerosol = ', iaero_n2 64 !!read altitudes and haze mmrs 65 integer Nfine 66 !parameter(Nfine=21) 67 parameter(Nfine=701) 68 character(len=100) :: file_path 69 real,save :: levdat(Nfine),densdat(Nfine) 82 70 83 if (dusttau.gt.0) then 84 ia=ia+1 85 iaero_dust=ia 86 endif 87 if (is_master) write(*,*) '--- Dust aerosol = ', iaero_dust 71 !---------------- INPUT ------------------------------------------------ 88 72 89 90 if (aeroback2lay) then 91 ia=ia+1 92 iaero_back2lay=ia 93 endif 94 if (is_master) write(*,*) '--- Two-layer aerosol = ', iaero_back2lay 73 !! Read data 74 IF (firstcall) then 75 firstcall=.false. 76 file_path=trim(datadir)//'/haze_prop/hazemmr.txt' 77 open(224,file=file_path,form='formatted') 78 do ifine=1,Nfine 79 read(224,*) levdat(ifine), densdat(ifine) 80 enddo 81 close(224) 82 print*, 'TB22 read Haze MMR profile' 83 ENDIF 95 84 96 if (aeronh3) then 97 ia=ia+1 98 iaero_nh3=ia 99 endif 100 if (is_master) write(*,*) '--- NH3 Cloud = ', iaero_nh3 85 !! Interpolate on the model vertical grid 86 do ig=1,ngrid 87 CALL interp_line(levdat,densdat,Nfine,zzlay(ig,:)/1000.,profmmr(ig,:),nlayer) 88 enddo 101 89 102 if (aeronlay) then 103 do i=1,nlayaero 104 ia=ia+1 105 iaero_nlay(i)=ia 106 enddo 107 endif 108 if (is_master) write(*,*) '--- N-layer aerosol = ', iaero_nlay 90 !! Get profile Mass mixing ratio from number density: part.cm-3 --> m-3 --> m3 m-3 91 ! --> kg m-3 --> kg/kg 92 do iaer=1,naerkind 93 if(iaer.eq.iaero_haze.and.1.eq.2) then !TB22 activate/deactivate mmr or part density 94 !print*, 'Haze profile is fixed' 95 do ig=1,ngrid 96 do l=1,nlayer 97 !from number density in cm-3 98 profmmr(ig,l)=profmmr(ig,l)*1.e6*4./3.*pi*reffrad(ig,l,iaer)**3*rho_q(i_haze)/(pplay(ig,l)/(r*pt(ig,l))) 99 enddo 100 enddo 101 endif 102 enddo 103 end subroutine haze_prof 104 !================================================================== 109 105 110 if (aeroaurora) then111 ia=ia+1112 iaero_aurora=ia113 endif114 if (is_master) write(*,*) '--- Auroral aerosols = ', iaero_aurora115 116 if (aerogeneric .ne. 0) then117 do i=1,aerogeneric118 ia = ia+1119 iaero_generic(i) = ia120 enddo121 endif122 123 if (is_master) then124 write(*,*)'--- Radiative Generic Condensable Species = ',iaero_generic125 126 write(*,*) '=== Number of aerosols= ', ia127 endif ! of is_master128 129 ! For the zero aerosol case, we currently make a dummy n2 aerosol which is zero everywhere.130 ! (See aeropacity.F90 for how this works). A better solution would be to turn off the131 ! aerosol machinery in the no aerosol case, but this would be complicated. LK132 133 if (ia.eq.0) then !For the zero aerosol case.134 ia = 0135 noaero = .true.136 endif137 138 if (.not.noaero .and. ia.ne.naerkind) then139 if (is_master) then140 print*, 'Aerosols counted not equal to naerkind'141 print*, 'set correct value for nearkind in callphys.def'142 print*, 'which should be ',ia143 print*, 'according to current options in callphys.def'144 print*, 'or change/correct incompatible options there'145 print*, 'Abort in iniaerosol'146 endif147 call abort_physic("iniaerosl",'wrong number of aerosols',1)148 endif ! of if (ia.ne.naerkind)149 150 END SUBROUTINE iniaerosol151 106 152 107 end module aerosol_mod -
trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90
r3193 r3195 6 6 7 7 subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, & 8 albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & 8 albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & 9 9 tsurf,fract,dist_star,aerosol,muvar, & 10 10 dtlw,dtsw,fluxsurf_lw, & … … 26 26 use ioipsl_getin_p_mod, only: getin_p 27 27 use gases_h, only: ngasmx 28 use radii_mod, only : su_aer_radii,n2_reffrad,dust_reffrad,h2so4_reffrad,back2lay_reffrad 29 use aerosol_mod, only : iaero_n2,iaero_dust,iaero_h2so4, & 30 iaero_back2lay, iaero_aurora 28 use radii_mod, only : su_aer_radii 29 use aerosol_mod, only : iaero_haze 31 30 use aeropacity_mod, only: aeropacity 32 31 use aeroptproperties_mod, only: aeroptproperties 33 use tracer_h, only: igcm_n2_ice34 32 use tracer_h, only: constants_epsi_generic 35 33 use comcstfi_mod, only: pi, mugaz, cpp … … 37 35 diagdtau,kastprof,strictboundcorrk,specOLR, & 38 36 tplanckmin,tplanckmax,global1d, & 39 generic_condensation 37 generic_condensation, aerohaze, haze_radproffix 40 38 use optcv_mod, only: optcv 41 39 use optci_mod, only: optci … … 57 55 ! 58 56 ! Authors 59 ! ------- 57 ! ------- 60 58 ! Emmanuel 01/2001, Forget 09/2001 61 59 ! Robin Wordsworth (2009) … … 65 63 !----------------------------------------------------------------------- 66 64 ! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid 67 ! Layer #1 is the layer near the ground. 65 ! Layer #1 is the layer near the ground. 68 66 ! Layer #nlayer is the layer at the top. 69 67 !----------------------------------------------------------------------- … … 90 88 logical,intent(in) :: firstcall ! Signals first call to physics. 91 89 logical,intent(in) :: lastcall ! Signals last call to physics. 92 90 93 91 ! OUTPUT 94 92 REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght. … … 109 107 REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! VI optical thickness of layers within narrowbands for diags (). 110 108 REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! IR optical thickness of layers within narrowbands for diags (). 111 112 113 114 115 116 ! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h. 109 110 111 112 113 114 ! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h. 117 115 ! made "save" variables so they are allocated once in for all, not because 118 116 ! the values need be saved from a time step to the next … … 208 206 real vtmp(nlayer) 209 207 REAL*8 muvarrad(L_LEVELS) 210 211 ! Included by MT for albedo calculations. 208 209 ! Included by MT for albedo calculations. 212 210 REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation. 213 211 REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation. 214 212 215 213 ! local variable 216 214 integer ok ! status (returned by NetCDF functions) … … 254 252 allocate(tauaero(L_LEVELS,naerkind)) 255 253 endif 256 257 if(.not.allocated(QXVAER)) then 254 255 if(.not.allocated(QXVAER)) then 258 256 allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind), stat=ok) 259 257 if (ok /= 0) then … … 334 332 335 333 call su_aer_radii(ngrid,nlayer,reffrad,nueffrad) 336 337 334 335 338 336 !-------------------------------------------------- 339 337 ! Set up correlated k … … 353 351 ! call sugas_corrk ! Set up gaseous absorption properties. 354 352 ! call suaer_corrk ! Set up aerosol optical properties. 355 353 356 354 357 355 ! now that L_NGAUSS has been initialized (by sugas_corrk) … … 435 433 endif 436 434 437 435 438 436 if(varfixed .and. generic_condensation)then 439 437 write(*,*) "Deep generic tracer vapor mixing ratio ? (no effect if negative) " … … 449 447 450 448 !======================================================================= 451 ! I.b Initialization on every call 449 ! I.b Initialization on every call 452 450 !======================================================================= 453 451 454 452 qxvaer(:,:,:)=0.0 455 453 qsvaer(:,:,:)=0.0 … … 467 465 ! Effective radius and variance of the aerosols 468 466 !-------------------------------------------------- 469 470 do iaer=1,naerkind 471 472 if ((iaer.eq.iaero_n2).and.tracer.and.(igcm_n2_ice.gt.0)) then ! Treat condensed n2 particles. 473 call n2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_n2)) 474 475 call planetwide_maxval(reffrad(:,:,iaero_n2),maxvalue) 476 call planetwide_minval(reffrad(:,:,iaero_n2),minvalue) 477 if (is_master) then 478 print*,'Max. N2 ice particle size = ',maxvalue/1.e-6,' um' 479 print*,'Min. N2 ice particle size = ',minvalue/1.e-6,' um' 480 end if 481 end if 482 ! Trear other aerosols here 483 484 if(iaer.eq.iaero_dust)then 485 call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust)) 486 if (is_master) then 487 print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um' 488 end if 489 endif 490 491 if(iaer.eq.iaero_h2so4)then 492 call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4)) 493 if (is_master) then 494 print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um' 495 end if 496 endif 497 498 if(iaer.eq.iaero_back2lay)then 499 call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev) 500 endif 501 502 ! For n-layer aerosol size set once for all at firstcall in su_aer_radii 503 504 ! if(iaer.eq.iaero_aurora)then 505 ! call aurora_reffrad(ngrid,nlayer,reffrad(1,1,iaero_aurora)) 506 ! endif 507 508 end do !iaer=1,naerkind. 467 if (aerohaze) then 468 do iaer=1,naerkind 469 if ((iaer.eq.iaero_haze)) then 470 call su_aer_radii(ngrid,nlayer,reffrad(1,1,iaer), & 471 nueffrad(1,1,iaer)) 472 endif 473 end do !iaer=1,naerkind. 474 if (haze_radproffix) then 475 print*, 'haze_radproffix=T : fixed profile for haze rad' 476 else 477 print*,'reffrad haze:',reffrad(1,1,iaero_haze) 478 print*,'nueff haze',nueffrad(1,1,iaero_haze) 479 endif 480 endif 509 481 510 482 … … 518 490 QVISsQREF3d,omegaVIS3d,gVIS3d, & 519 491 QIRsQREF3d,omegaIR3d,gIR3d, & 520 QREFvis3d,QREFir3d) 492 QREFvis3d,QREFir3d) 521 493 522 494 ! Get aerosol optical depths. 523 495 call aeropacity(ngrid,nlayer,nq,pplay,pplev, pt,pq,aerosol, & 524 reffrad,nueffrad,QREFvis3d,QREFir3d, & 525 tau_col,cloudfrac,totcloudfrac,clearsky) 526 527 !----------------------------------------------------------------------- 496 reffrad,nueffrad,QREFvis3d,QREFir3d, & 497 tau_col,cloudfrac,totcloudfrac,clearsky) 498 499 !----------------------------------------------------------------------- 528 500 do ig=1,ngrid ! Starting Big Loop over every GCM column 529 501 !----------------------------------------------------------------------- … … 539 511 ! The transformation in the vertical is the same as for temperature. 540 512 !----------------------------------------------------------------------- 541 542 513 514 543 515 do iaer=1,naerkind 544 516 ! Shortwave. 545 do nw=1,L_NSPECTV 546 517 do nw=1,L_NSPECTV 518 547 519 do l=1,nlayer 548 520 … … 580 552 581 553 end do ! L_NSPECTV 582 554 583 555 do nw=1,L_NSPECTI 584 556 ! Longwave … … 618 590 619 591 end do ! L_NSPECTI 620 592 621 593 end do ! naerkind 622 594 … … 627 599 do nw=1,L_NSPECTV 628 600 if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then 629 message='Serious problems with qsvaer values' 601 message='Serious problems with qsvaer values' 630 602 call abort_physic(subname,message,1) 631 603 endif … … 635 607 end do 636 608 637 do nw=1,L_NSPECTI 609 do nw=1,L_NSPECTI 638 610 if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then 639 message='Serious problems with qsvaer values' 611 message='Serious problems with qsvaer values' 640 612 call abort_physic(subname,message,1) 641 613 endif … … 651 623 ! Aerosol optical depths 652 624 !----------------------------------------------------------------------- 653 654 do iaer=1,naerkind ! a bug was here 625 626 do iaer=1,naerkind ! a bug was here 655 627 do k=0,nlayer-1 656 628 657 629 pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ & 658 630 (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1)) … … 667 639 tauaero(1,iaer) = tauaero(2,iaer) 668 640 !tauaero(1,iaer) = 0. 669 !JL18 at time of testing, the two above conditions gave the same results bit for bit. 670 641 !JL18 at time of testing, the two above conditions gave the same results bit for bit. 642 671 643 end do ! naerkind 672 644 … … 696 668 697 669 call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic) 698 670 699 671 if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer 700 672 … … 704 676 do l=1,nlayer 705 677 qvar(2*l) = pq(ig,nlayer+1-l,i_var) 706 qvar(2*l+1) = pq(ig,nlayer+1-l,i_var) 707 !JL13index qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2 678 qvar(2*l+1) = pq(ig,nlayer+1-l,i_var) 679 !JL13index qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2 708 680 !JL13index ! Average approximation as for temperature... 709 681 end do … … 713 685 714 686 do l=1,nlayer ! Here we will assign fixed water vapour profiles globally. 715 687 716 688 call Psat_generic(pt(ig,l),pplay(ig,l),metallicity,psat,qsat) 717 689 718 690 if (qsat .lt. qvap_deep) then 719 691 pq_temp(l) = qsat ! fully saturated everywhere 720 else 692 else 721 693 pq_temp(l) = qvap_deep 722 694 end if 723 695 724 696 end do 725 697 726 698 do l=1,nlayer 727 699 qvar(2*l) = pq_temp(nlayer+1-l) 728 700 qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2 729 701 end do 730 702 731 703 qvar(1)=qvar(2) 732 704 733 705 else 734 706 do k=1,L_LEVELS … … 739 711 endif 740 712 741 end do ! do iq=1,nq loop on tracers 713 end do ! do iq=1,nq loop on tracers 742 714 743 715 end if ! if (generic_condensation) … … 780 752 muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2 781 753 END DO 782 754 783 755 muvarrad(1) = muvarrad(2) 784 muvarrad(2*nlayer+1)=muvar(ig,1) 785 756 muvarrad(2*nlayer+1)=muvar(ig,1) 757 786 758 ! Keep values inside limits for which we have radiative transfer coefficients !!! 787 759 if(L_REFVAR.gt.1)then ! (there was a bug here) … … 805 777 tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2 806 778 END DO 807 779 808 780 plevrad(1) = 0. 809 ! plevrad(2) = 0. !! JL18 enabling this line puts the radiative top at p=0 which was the idea before, but does not seem to perform best after all. 810 781 ! plevrad(2) = 0. !! JL18 enabling this line puts the radiative top at p=0 which was the idea before, but does not seem to perform best after all. 782 811 783 tlevrad(1) = tlevrad(2) 812 784 tlevrad(2*nlayer+1)=tsurf(ig) 813 814 pmid(1) = pplay(ig,nlayer)/scalep 785 786 pmid(1) = pplay(ig,nlayer)/scalep 815 787 pmid(2) = pmid(1) 816 788 817 789 tmid(1) = tlevrad(2) 818 790 tmid(2) = tmid(1) 819 791 820 792 DO l=1,L_NLAYRAD-1 821 793 tmid(2*l+1) = tlevrad(2*l+1) … … 829 801 !!Alternative interpolation: 830 802 ! pmid(3) = pmid(1) 831 ! pmid(4) = pmid(1) 803 ! pmid(4) = pmid(1) 832 804 ! tmid(3) = tmid(1) 833 805 ! tmid(4) = tmid(1) … … 870 842 else 871 843 print*,'***********************************************' 872 print*,'we allow model to continue with tlevrad<tgasmin' 844 print*,'we allow model to continue with tlevrad<tgasmin' 873 845 print*,' ... we assume we know what you are doing ... ' 874 846 print*,' ... but do not let this happen too often ... ' … … 886 858 else 887 859 print*,'***********************************************' 888 print*,'we allow model to continue with tlevrad>tgasmax' 860 print*,'we allow model to continue with tlevrad>tgasmax' 889 861 print*,' ... we assume we know what you are doing ... ' 890 862 print*,' ... but do not let this happen too often ... ' … … 1000 972 1001 973 ! Equivalent Albedo Calculation (for OUTPUT). MT2015 1002 if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight. 1003 surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV)) 974 if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight. 975 surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV)) 1004 976 if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive. 1005 DO nw=1,L_NSPECTV 977 DO nw=1,L_NSPECTV 1006 978 albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw) 1007 979 ENDDO … … 1025 997 1026 998 call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, & 1027 wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, & 999 wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, & 1028 1000 fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi) 1029 1001 … … 1033 1005 1034 1006 ! Flux incident at the top of the atmosphere 1035 fluxtop_dn(ig)=fluxtopvdn 1007 fluxtop_dn(ig)=fluxtopvdn 1036 1008 1037 1009 fluxtop_lw(ig) = real(nfluxtopi) … … 1039 1011 fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD)) 1040 1012 fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD)) 1041 1042 ! Flux absorbed by the surface. By MT2015. 1013 1014 ! Flux absorbed by the surface. By MT2015. 1043 1015 fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig)) 1044 1016 … … 1056 1028 ! Spectral output, for exoplanet observational comparison 1057 1029 if(specOLR)then 1058 do nw=1,L_NSPECTI 1030 do nw=1,L_NSPECTI 1059 1031 OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth 1060 1032 end do 1061 do nw=1,L_NSPECTV 1033 do nw=1,L_NSPECTV 1062 1034 GSR_nu(ig,nw)=nfluxgndv_nu(nw)/DWNV(nw) 1063 1035 OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth … … 1072 1044 dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) & 1073 1045 *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) 1074 END DO 1046 END DO 1075 1047 1076 1048 ! These are values at top of atmosphere … … 1098 1070 enddo 1099 1071 enddo 1100 endif 1101 1102 1103 !----------------------------------------------------------------------- 1072 endif 1073 1074 1075 !----------------------------------------------------------------------- 1104 1076 end do ! End of big loop over every GCM column. 1105 1077 !----------------------------------------------------------------------- … … 1118 1090 open(116,file='surf_vals.out') 1119 1091 write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), & 1120 real(-nfluxtopv),real(nfluxtopi) 1092 real(-nfluxtopv),real(nfluxtopi) 1121 1093 close(116) 1122 1094 … … 1147 1119 write(118,*) plevrad(2*l) 1148 1120 do nw=1,L_NSPECTI 1149 write(119,*) fluxupi_nu(l,nw) 1121 write(119,*) fluxupi_nu(l,nw) 1150 1122 enddo 1151 enddo 1123 enddo 1152 1124 close(118) 1153 1125 close(119) -
trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90
r3193 r3195 1 1 MODULE callkeys_mod 2 IMPLICIT NONE 2 IMPLICIT NONE 3 3 4 4 logical,save :: callrad,corrk,calldifv,UseTurbDiff … … 10 10 logical,save :: callgasvis,continuum,graybody 11 11 !$OMP THREADPRIVATE(callgasvis,continuum,graybody) 12 logical,save :: strictboundcorrk 12 logical,save :: strictboundcorrk 13 13 !$OMP THREADPRIVATE(strictboundcorrk) 14 logical,save :: strictboundcia 14 logical,save :: strictboundcia 15 15 !$OMP THREADPRIVATE(strictboundcia) 16 16 … … 57 57 58 58 !! Pluto-specific variables 59 logical,save :: haze_radproffix,haze_proffix 60 !$OMP THREADPRIVATE(haze_radproffix,haze_proffix) 61 logical,save :: haze,fasthaze,changeti,changetid,aerohaze 62 !$OMP THREADPRIVATE(haze,fasthaze,changeti,changetid,aerohaze) 63 logical,save :: fast,metcloud,monoxcloud,glaflow,triton,paleo,nlte,strobel 64 !$OMP THREADPRIVATE(fast,metcloud,monoxcloud,glaflow,triton,paleo,nlte,strobel 65 logical,save :: kbo,nbsub,cooling,mode_n2,thres_n2ice,thres_coice 66 !$OMP THREADPRIVATE(kbo,nbsub,cooling,mode_n2,thres_n2ice,thres_coice) 67 logical,save :: deltab,nb_monomer 68 !$OMP THREADPRIVATE(deltab,nb_monomer) 69 logical,save :: metlateq,tholateq,metls1,metls2 70 !$OMP THREADPRIVATE(metlateq,tholateq,metls1,metls2) 71 logical,save :: mode_ch4,mode_tholins 72 !$OMP THREADPRIVATE(mode_ch4,mode_tholins) 73 logical,save :: source_haze,ch4lag,latlag,vmrlag,kfix,hazeconservch4 74 !$OMP THREADPRIVATE(source_haze,ch4lag,latlag,vmrlag,kfix,hazeconservch4) 75 logical,save :: fracsource,latsource,lonsource 76 !$OMP THREADPRIVATE(fracsource,latsource,lonsource) 77 logical,save :: spelon1,spelon2,spelat1,spelat2,specalb 78 !$OMP THREADPRIVATE(spelon1,spelon2,spelat1,spelat2,specalb) 79 logical,save :: assymflux,mode_hs,tsurfmax,albmin_ch4 80 !$OMP THREADPRIVATE(assymflux,mode_hs,tsurfmax,albmin_ch4) 81 logical,save :: feedback_met,thres_ch4ice,fdch4_latn,fdch4_lats 82 !$OMP THREADPRIVATE(feedback_met,thres_ch4ice,fdch4_latn,fdch4_lats) 83 logical,save :: globmean1d,kmix_proffix,rad_haze 84 !$OMP THREADPRIVATE(globmean1d,kmix_proffix,rad_haze) 85 logical,save :: fdch4_lone,fdch4_lonw,fdch4_maxalb,fdch4_depalb,fdch4_finalb 86 !$OMP THREADPRIVATE(fdch4_lone,fdch4_lonw,fdch4_maxalb,fdch4_depalb,fdch4_finalb 87 logical,save :: tholatn,tholats,tholone,tholonw 88 !$OMP THREADPRIVATE(tholatn,tholats,tholone,tholonw) 89 logical,save :: fdch4_ampl,fdch4_maxice,condmetsurf,condcosurf,vertdiff 90 !$OMP THREADPRIVATE(fdch4_ampl,fdch4_maxice,condmetsurf,condcosurf,vertdiff) 91 logical,save :: convergeps,conservn2,patchflux,condensn2,no_n2frost 92 !$OMP THREADPRIVATE(convergeps,conservn2,patchflux,condensn2,no_n2frost 93 59 logical,save :: methane,carbox 60 !$OMP THREADPRIVATE(methane,carbox) 61 logical,save :: haze,haze_proffix,haze_radproffix 62 !$OMP THREADPRIVATE(haze,haze_proffix,haze_radproffix) 63 logical,save :: fasthaze,changeti,changetid,aerohaze,fractal 64 !$OMP THREADPRIVATE(fasthaze,changeti,changetid,aerohaze,fractal) 65 logical,save :: fast,metcloud,monoxcloud,glaflow,triton,paleo 66 !$OMP THREADPRIVATE(fast,metcloud,monoxcloud,glaflow,triton,paleo) 67 logical,save :: nlte,strobel 68 !$OMP THREADPRIVATE(nlte,strobel) 69 logical,save :: kbo 70 !$OMP THREADPRIVATE(kbo) 71 logical,save :: cooling 72 !$OMP THREADPRIVATE(cooling) 73 logical,save :: source_haze,hazeconservch4 74 !$OMP THREADPRIVATE(source_haze,hazeconservch4) 75 logical,save :: ch4lag,tsurfmax 76 !$OMP THREADPRIVATE(ch4lag,tsurfmax) 77 logical,save :: specalb 78 !$OMP THREADPRIVATE(specalb) 79 logical,save :: assymflux 80 !$OMP THREADPRIVATE(assymflux) 81 logical,save :: condmetsurf,condcosurf,vertdiff 82 !$OMP THREADPRIVATE(condmetsurf,condcosurf,vertdiff) 83 logical,save :: convergeps,conservn2,condensn2,no_n2frost 84 !$OMP THREADPRIVATE(convergeps,conservn2,condensn2,no_n2frost) 85 logical,save :: globmean1d,kmix_proffix 86 !$OMP THREADPRIVATE(globmean1d,kmix_proffix) 87 integer,save :: nbsub 88 !$OMP THREADPRIVATE(nbsub) 89 integer,save :: mode_n2 90 !$OMP THREADPRIVATE(mode_n2) 91 integer,save :: mode_ch4 92 !$OMP THREADPRIVATE(mode_ch4) 93 integer,save :: mode_tholins 94 !$OMP THREADPRIVATE(mode_tholins) 95 integer,save :: kfix 96 !$OMP THREADPRIVATE(kfix) 97 integer,save :: mode_hs 98 !$OMP THREADPRIVATE(mode_hs) 99 integer,save :: feedback_met 100 !$OMP THREADPRIVATE(feedback_met) 101 integer,save :: patchflux 102 !$OMP THREADPRIVATE(patchflux) 103 integer,save :: nb_monomer 104 !$OMP THREADPRIVATE(nb_monomer) 105 real,save :: Nmix_co 106 !$OMP THREADPRIVATE(Nmix_co) 107 real,save :: Nmix_ch4 108 !$OMP THREADPRIVATE(Nmix_ch4) 109 real,save :: tau_n2 110 !$OMP THREADPRIVATE(tau_n2) 111 real,save :: tau_ch4 112 !$OMP THREADPRIVATE(tau_ch4) 113 real,save :: tau_co 114 !$OMP THREADPRIVATE(tau_co) 115 real,save :: tau_prechaze 116 !$OMP THREADPRIVATE(tau_prechaze) 117 real,save :: paleoyears 118 !$OMP THREADPRIVATE(paleoyears) 119 real,save :: dayfrac 120 !$OMP THREADPRIVATE(dayfrac) 121 real,save :: thresh_non2 122 !$OMP THREADPRIVATE(thresh_non2) 123 real,save :: vmrch4fix 124 !$OMP THREADPRIVATE(vmrch4fix) 125 real,save :: fluxgeo 126 !$OMP THREADPRIVATE(fluxgeo) 127 real,save :: fluxgeo2 128 !$OMP THREADPRIVATE(fluxgeo2) 129 real,save :: deltab 130 !$OMP THREADPRIVATE(deltab) 131 real,save :: metlateq 132 !$OMP THREADPRIVATE(metlateq) 133 real,save :: metls1 134 !$OMP THREADPRIVATE(metls1) 135 real,save :: metls2 136 !$OMP THREADPRIVATE(metls2) 137 real,save :: tholateq 138 !$OMP THREADPRIVATE(tholateq) 139 real,save :: tholatn 140 !$OMP THREADPRIVATE(tholatn) 141 real,save :: tholats 142 !$OMP THREADPRIVATE(tholats) 143 real,save :: tholone 144 !$OMP THREADPRIVATE(tholone) 145 real,save :: tholonw 146 !$OMP THREADPRIVATE(tholonw) 147 real,save :: spelon1 148 !$OMP THREADPRIVATE(spelon1) 149 real,save :: spelat2 150 !$OMP THREADPRIVATE(spelat2) 151 real,save :: spelat1 152 !$OMP THREADPRIVATE(spelat1) 153 real,save :: spelon2 154 !$OMP THREADPRIVATE(spelon2) 155 real,save :: latlag 156 !$OMP THREADPRIVATE(latlag) 157 real,save :: vmrlag 158 !$OMP THREADPRIVATE(vmrlag) 159 real,save :: albmin_ch4 160 !$OMP THREADPRIVATE(albmin_ch4) 161 real,save :: fracsource 162 !$OMP THREADPRIVATE(fracsource) 163 real,save :: latsource 164 !$OMP THREADPRIVATE(latsource) 165 real,save :: lonsource 166 !$OMP THREADPRIVATE(lonsource) 167 real,save :: thres_ch4ice 168 !$OMP THREADPRIVATE(thres_ch4ice) 169 real,save :: thres_n2ice 170 !$OMP THREADPRIVATE(thres_n2ice) 171 real,save :: thres_coice 172 !$OMP THREADPRIVATE(thres_coice) 173 real,save :: fdch4_latn 174 !$OMP THREADPRIVATE(fdch4_latn) 175 real,save :: fdch4_lats 176 !$OMP THREADPRIVATE(fdch4_lats) 177 real,save :: fdch4_lone 178 !$OMP THREADPRIVATE(fdch4_lone) 179 real,save :: fdch4_lonw 180 !$OMP THREADPRIVATE(fdch4_lonw) 181 real,save :: fdch4_maxice 182 !$OMP THREADPRIVATE(fdch4_maxice) 183 real,save :: fdch4_maxalb 184 !$OMP THREADPRIVATE(fdch4_maxalb) 185 real,save :: fdch4_ampl 186 !$OMP THREADPRIVATE(fdch4_ampl) 187 real,save :: fdch4_depalb 188 !$OMP THREADPRIVATE(fdch4_depalb) 189 real,save :: fdch4_finalb 190 !$OMP THREADPRIVATE(fdch4_finalb) 191 real,save :: rad_haze 192 !$OMP THREADPRIVATE(rad_haze) 94 193 95 194 logical,save :: global1d … … 169 268 real,save :: kmixmin 170 269 !$OMP THREADPRIVATE(kmixmin) 171 270 172 271 logical,save :: iscallphys=.false.!existence of callphys.def 173 272 !$OMP THREADPRIVATE(iscallphys) 174 273 175 274 ! do we read a startphy.nc file (default=.true.) 176 logical,save :: startphy_file=.true. 275 logical,save :: startphy_file=.true. 177 276 !$OMP THREADPRIVATE(startphy_file) 178 277 -
trunk/LMDZ.PLUTO/libf/phypluto/callsedim.F
r3193 r3195 6 6 ! use radii_mod, only: h2o_reffrad 7 7 ! use aerosol_mod, only : iaero_h2o 8 USE tracer_h, only : igcm_ n2_ice,radius,rho_q8 USE tracer_h, only : igcm_haze,radius,rho_q 9 9 use comcstfi_mod, only: g 10 10 … … 12 12 13 13 !================================================================== 14 ! 14 ! 15 15 ! Purpose 16 16 ! ------- 17 17 ! Calculates sedimentation of aerosols depending on their 18 18 ! density and radius. 19 ! 19 ! 20 20 ! Authors 21 21 ! ------- 22 22 ! F. Forget (1999) 23 23 ! Tracer generalisation by E. Millour (2009) 24 ! 24 ! 25 25 !================================================================== 26 26 … … 43 43 real,intent(in) :: pdqfi(ngrid,nlay,nq) ! tendency on tracers before 44 44 ! sedimentation (kg/kg.s-1) 45 45 46 46 real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1) 47 47 real,intent(out) :: pdqs_sed(ngrid,nq) ! flux at surface (kg.m-2.s-1) … … 53 53 54 54 ! for particles with varying radii: 55 real,allocatable,save :: reffrad(:,:,:) ! particle radius (m) 55 real,allocatable,save :: reffrad(:,:,:) ! particle radius (m) 56 56 real,allocatable,save :: nueffrad(:,:,:) ! aerosol effective radius variance 57 57 !$OMP THREADPRIVATE(reffrad,nueffrad) … … 80 80 allocate(nueffrad(ngrid,nlay,naerkind)) 81 81 ENDIF ! of IF (firstcall) 82 82 83 83 !======================================================================= 84 84 ! Preliminary calculation of the layer characteristics … … 87 87 do l=1,nlay 88 88 do ig=1, ngrid 89 masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 89 masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 90 90 epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l) 91 91 zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep … … 97 97 ! [This has been rearranged by L. Kerber to allow the sedimentation 98 98 ! of general tracers.] 99 99 100 100 do iq=1,nq 101 if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_ n2_ice)) then ! JVO 08/2017 : be careful radius was tested uninitialized (fixed) ...102 103 ! (no sedim for gases, and n2_ice sedim is done in condense_n2) 101 if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_haze)) then ! JVO 08/2017 : be careful radius was tested uninitialized (fixed) ... 102 103 ! (no sedim for gases, and n2_ice sedim is done in condense_n2) 104 104 105 105 ! store locally updated tracers 106 106 107 do l=1,nlay 107 do l=1,nlay 108 108 do ig=1, ngrid 109 109 zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep 110 110 enddo 111 111 enddo ! of do l=1,nlay 112 112 113 113 !====================================================================== 114 ! Sedimentation 114 ! Sedimentation 115 115 !====================================================================== 116 116 ! Water !AF24: deleted 117 117 118 118 ! ! General Case 119 ! else 119 ! else 120 120 call newsedim(ngrid,nlay,1,ptimestep, 121 121 & pplev,masse,epaisseur,zt,radius(iq),rho_q(iq), -
trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90
r3193 r3195 1 subroutine condens _n2(klon,nlayer,nq,ptimestep, &1 subroutine condense_n2(klon,klev,nq,ptimestep, & 2 2 pcapcal,pplay,pplev,ptsrf,pt, & 3 3 pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, & … … 8 8 use radinc_h, only : naerkind 9 9 use comgeomfi_h 10 use comcstfi_mod, only: g, r, cpp 11 USE surfdat_h, only: phisfi 12 USE tracer_h, only: noms, igcm_n2 13 USE callkeys_mod, only: fast,ch4lag,latlag,lw_n2,nbsub,no_n2frost,tsurfmax,kmixmin,source_haze,vmrlag 14 USE dimphy, only : klon,klev 10 use comcstfi_mod, only: g, r, cpp, pi 11 USE surfdat_h, only: phisfi,kp,p00 12 USE tracer_h, only: noms, igcm_n2, lw_n2 13 USE callkeys_mod, only: fast,ch4lag,latlag,nbsub,no_n2frost,tsurfmax,kmixmin,source_haze,vmrlag 14 USE comvert_mod, ONLY: ap,bp 15 use geometry_mod, only: latitude 15 16 16 17 … … 22 23 ! Condense and/or sublime N2 ice on the ground and in the 23 24 ! atmosphere, and sediment the ice. 24 ! 25 ! 25 26 ! Inputs 26 ! ------ 27 ! ------ 27 28 ! klon Number of vertical columns 28 ! nlayerNumber of layers29 ! pplay(klon, nlayer) Pressure layers30 ! pplev(klon, nlayer+1) Pressure levels31 ! pt(klon, nlayer) Temperature (in K)29 ! klev Number of layers 30 ! pplay(klon,klev) Pressure layers 31 ! pplev(klon,klev+1) Pressure levels 32 ! pt(klon,klev) Temperature (in K) 32 33 ! ptsrf(klon) Surface temperature 33 ! 34 ! pdt(klon, nlayermx) Time derivative before condensation/sublimation of pt34 ! 35 ! pdt(klon,klev) Time derivative before condensation/sublimation of pt 35 36 ! pdtsrf(klon) Time derivative before condensation/sublimation of ptsrf 36 37 ! picen2(klon) n2 ice at the surface (kg/m2) 37 ! 38 ! 38 39 ! Outputs 39 40 ! ------- 40 41 ! pdpsrf(klon) \ Contribution of condensation/sublimation 41 ! pdtc(klon, nlayermx) / to the time derivatives of Ps, pt, and ptsrf42 ! pdtc(klon,klev) / to the time derivatives of Ps, pt, and ptsrf 42 43 ! pdtsrfc(klon) / 43 44 ! pdicen2(klon) Tendancy n2 ice at the surface (kg/m2) 44 ! 45 ! 45 46 ! Both 46 47 ! ---- … … 49 50 ! 50 51 ! Authors 51 ! ------- 52 ! ------- 52 53 ! Francois Forget (1996,2013) 53 54 ! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009) 54 ! 55 ! 55 ! 56 ! 56 57 !================================================================== 57 58 … … 59 60 ! Arguments 60 61 61 INTEGER klon, nlayer, nq62 63 REAL ptimestep 64 REAL pplay(klon, nlayer),pplev(klon,nlayer+1)62 INTEGER klon, klev, nq 63 64 REAL ptimestep 65 REAL pplay(klon,klev),pplev(klon,klev+1) 65 66 REAL pcapcal(klon) 66 REAL pt(klon, nlayer)67 REAL pt(klon,klev) 67 68 REAL ptsrf(klon),flu1(klon),flu2(klon),flu3(klon) 68 REAL pphi(klon, nlayer)69 REAL pdt(klon, nlayer),pdtsrf(klon),pdtc(klon,nlayer)69 REAL pphi(klon,klev) 70 REAL pdt(klon,klev),pdtsrf(klon),pdtc(klon,klev) 70 71 REAL pdtsrfc(klon),pdpsrf(klon) 71 72 REAL picen2(klon),psolaralb(klon),pemisurf(klon) 72 73 74 REAL pu(klon, nlayer) , pv(klon,nlayer)75 REAL pdu(klon, nlayer) , pdv(klon,nlayer)76 REAL pduc(klon, nlayer) , pdvc(klon,nlayer)77 REAL pq( ngridmx,nlayer,nq),pdq(klon,nlayer,nq)78 REAL pdqc(klon, nlayer,nq)73 74 75 REAL pu(klon,klev) , pv(klon,klev) 76 REAL pdu(klon,klev) , pdv(klon,klev) 77 REAL pduc(klon,klev) , pdvc(klon,klev) 78 REAL pq(klon,klev,nq),pdq(klon,klev,nq) 79 REAL pdqc(klon,klev,nq) 79 80 80 81 !----------------------------------------------------------------------- … … 83 84 INTEGER l,ig,ilay,it,iq,ich4_gas 84 85 85 REAL*8 zt( ngridmx,nlayermx)86 REAL*8 zt(klon,klev) 86 87 REAL tcond_n2 87 88 REAL pcond_n2 88 89 REAL glob_average2d ! function 89 REAL zqn2( ngridmx,nlayermx) ! N2 MMR used to compute Tcond/zqn290 REAL ztcond ( ngridmx,nlayermx)91 REAL ztcondsol( ngridmx),zfallheat92 REAL pdicen2( ngridmx)93 REAL zcondicea( ngridmx,nlayermx), zcondices(ngridmx)94 REAL zfallice( ngridmx,nlayermx+1)95 REAL zmflux( nlayermx+1)96 REAL zu( nlayermx),zv(nlayermx)97 REAL zq( nlayermx,nqmx),zq1(nlayermx)98 REAL ztsrf( ngridmx)99 REAL ztc( nlayermx), ztm(nlayermx+1)100 REAL zum( nlayermx+1) , zvm(nlayermx+1)101 REAL zqm( nlayermx+1,nqmx),zqm1(nlayermx+1)102 LOGICAL condsub( ngridmx)90 REAL zqn2(klon,klev) ! N2 MMR used to compute Tcond/zqn2 91 REAL ztcond (klon,klev) 92 REAL ztcondsol(klon),zfallheat 93 REAL pdicen2(klon) 94 REAL zcondicea(klon,klev), zcondices(klon) 95 REAL zfallice(klon,klev+1) 96 REAL zmflux(klev+1) 97 REAL zu(klev),zv(klev) 98 REAL zq(klev,nq),zq1(klev) 99 REAL ztsrf(klon) 100 REAL ztc(klev), ztm(klev+1) 101 REAL zum(klev+1) , zvm(klev+1) 102 REAL zqm(klev+1,nq),zqm1(klev+1) 103 LOGICAL condsub(klon) 103 104 REAL subptimestep 104 105 Integer Ntime 105 real masse ( nlayermx),w(nlayermx+1)106 real wq( ngridmx,nlayermx+1)106 real masse (klev),w(klev+1) 107 real wq(klon,klev+1) 107 108 real vstokes,reff 108 109 real dWtotsn2 109 real condnconsn2( ngridmx)110 real condnconsn2(klon) 110 111 real nconsMAXn2 111 112 ! Special diagnostic variables 112 real tconda1( ngridmx,nlayermx)113 real tconda2( ngridmx,nlayermx)114 real zdtsig ( ngridmx,nlayermx)115 real zdtlatent ( ngridmx,nlayermx)116 real zdt ( ngridmx,nlayermx)117 REAL albediceF( ngridmx)118 SAVE albediceF113 real tconda1(klon,klev) 114 real tconda2(klon,klev) 115 real zdtsig (klon,klev) 116 real zdtlatent (klon,klev) 117 real zdt (klon,klev) 118 REAL albediceF(klon) 119 ! SAVE albediceF 119 120 INTEGER nsubtimestep,itsub !number of subtimestep when calling vl1d 120 121 REAL subtimestep !ptimestep/nsubtimestep 121 REAL dtmax 122 123 REAL zplevhist( ngridmx)124 REAL zplevnew( ngridmx)125 REAL zplev( ngridmx)126 REAL zpicen2( ngridmx)127 REAL ztsrfhist( ngridmx)128 REAL zdtsrf( ngridmx)122 REAL dtmax 123 124 REAL zplevhist(klon) 125 REAL zplevnew(klon) 126 REAL zplev(klon) 127 REAL zpicen2(klon) 128 REAL ztsrfhist(klon) 129 REAL zdtsrf(klon) 129 130 REAL globzplevnew 130 131 131 REAL vmrn2( ngridmx)132 SAVE vmrn2133 REAL stephan 132 REAL vmrn2(klon) 133 ! SAVE vmrn2 134 REAL stephan 134 135 DATA stephan/5.67e-08/ ! Stephan Boltzman constant 135 136 SAVE stephan … … 167 168 IF (firstcall) THEN 168 169 ccond=cpp/(g*lw_n2) 169 print*,'In condens _n2cloud: ccond=',ccond,' latcond=',lw_n2170 print*,'In condense_n2cloud: ccond=',ccond,' latcond=',lw_n2 170 171 171 172 ! calculate global mean surface pressure for the fast mode 172 173 IF (fast) THEN 173 DO ig=1, ngridmx174 DO ig=1,klon 174 175 kp(ig) = exp(-phisfi(ig)/(r*38.)) 175 176 ENDDO … … 179 180 vmrn2(:) = 1. 180 181 IF (ch4lag) then 181 DO ig=1, ngridmx182 if (lati (ig)*180./pi.ge.latlag) then182 DO ig=1,klon 183 if (latitude(ig)*180./pi.ge.latlag) then 183 184 vmrn2(ig) = vmrlag 184 185 endif 185 186 ENDDO 186 ENDIF 187 ENDIF 187 188 IF (no_n2frost) then 188 DO ig=1, ngridmx189 DO ig=1,klon 189 190 if (picen2(ig).eq.0.) then 190 191 vmrn2(ig) = 1.e-15 191 192 endif 192 193 ENDDO 193 ENDIF 194 ENDIF 194 195 firstcall=.false. 195 196 ENDIF 196 197 197 198 !----------------------------------------------------------------------- 198 ! Calculation of n2 condensation / sublimation 199 ! Calculation of n2 condensation / sublimation 199 200 200 201 ! Variables used: 201 202 ! picen2(klon) : amount of n2 ice on the ground (kg/m2) 202 ! zcondicea(klon, nlayermx): condensation rate in layer l (kg/m2/s)203 ! zcondicea(klon,klev): condensation rate in layer l (kg/m2/s) 203 204 ! zcondices(klon) : condensation rate on the ground (kg/m2/s) 204 ! zfallice(klon, nlayermx) : amount of ice falling from layer l (kg/m2/s)205 ! zdtlatent(klon, nlayermx): dT/dt due to phase changes (K/s)205 ! zfallice(klon,klev) : amount of ice falling from layer l (kg/m2/s) 206 ! zdtlatent(klon,klev): dT/dt due to phase changes (K/s) 206 207 207 208 ! Tendencies initially set to 0 … … 213 214 pdicen2(1:klon) = 0. 214 215 zfallheat=0 215 pdqc(1:klon,1: nlayer,1:nq)=0216 pdtc(1:klon,1: nlayer)=0217 pduc(1:klon,1: nlayer)=0218 pdvc(1:klon,1: nlayer)=0219 zfallice(1:klon,1: nlayer+1)=0220 zcondicea(1:klon,1: nlayer)=0221 zdtlatent(1:klon,1: nlayer)=0222 zt(1:klon,1: nlayer)=0.216 pdqc(1:klon,1:klev,1:nq)=0 217 pdtc(1:klon,1:klev)=0 218 pduc(1:klon,1:klev)=0 219 pdvc(1:klon,1:klev)=0 220 zfallice(1:klon,1:klev+1)=0 221 zcondicea(1:klon,1:klev)=0 222 zdtlatent(1:klon,1:klev)=0 223 zt(1:klon,1:klev)=0. 223 224 224 225 !----------------------------------------------------------------------- … … 229 230 ! (calcul of zcondicea , zfallice and pdtc) 230 231 231 zt(1:klon,1: nlayer)=pt(1:klon,1:nlayer)+ pdt(1:klon,1:nlayer)*ptimestep232 zt(1:klon,1:klev)=pt(1:klon,1:klev)+ pdt(1:klon,1:klev)*ptimestep 232 233 if (igcm_n2.ne.0) then 233 zqn2(1:klon,1: nlayer) = 1. ! & temporaire234 ! zqn2(1:klon,1: nlayer)= pq(1:klon,1:nlayer,igcm_n2) + pdq(1:klon,1:nlayer,igcm_n2)*ptimestep234 zqn2(1:klon,1:klev) = 1. ! & temporaire 235 ! zqn2(1:klon,1:klev)= pq(1:klon,1:klev,igcm_n2) + pdq(1:klon,1:klev,igcm_n2)*ptimestep 235 236 else 236 zqn2(1:klon,1: nlayer) = 1.237 zqn2(1:klon,1:klev) = 1. 237 238 end if 238 239 239 240 if (.not.fast) then 240 241 ! Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2 241 DO l=1, nlayer242 DO l=1,klev 242 243 DO ig=1,klon 243 ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l)) 244 ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l)) 244 245 ENDDO 245 246 ENDDO 246 247 247 DO l= nlayer,1,-1248 DO l=klev,1,-1 248 249 DO ig=1,klon 249 250 pdtc(ig,l)=0. ! final tendancy temperature set to 0 … … 251 252 IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN 252 253 condsub(ig)=.true. !condensation at level l 253 IF (zfallice(ig,l+1).gt.0) then 254 zfallheat=zfallice(ig,l+1)*& 254 IF (zfallice(ig,l+1).gt.0) then 255 zfallheat=zfallice(ig,l+1)*& 255 256 (pphi(ig,l+1)-pphi(ig,l) +& 256 257 cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/lw_n2 … … 266 267 .AND. (zfallice(ig,l+1).gt.0)) THEN 267 268 268 zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/& 269 zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/& 269 270 (ccond*(pplev(ig,l)-pplev(ig,l+1))) 270 271 zcondicea(ig,l)= -zfallice(ig,l+1) … … 272 273 273 274 zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1) 274 275 275 276 END IF 276 277 277 278 ENDDO 278 279 ENDDO … … 285 286 ! Adding subtimesteps : in fast version, pressures too low lead to 286 287 ! instabilities. 287 IF (fast) THEN 288 IF (fast) THEN 288 289 IF (pplev(1,1).gt.0.3) THEN 289 nsubtimestep= 1 290 nsubtimestep= 1 290 291 ELSE 291 nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1) 292 nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1) 292 293 ENDIF 293 294 ELSE … … 307 308 ENDDO 308 309 ELSE 309 ! additional loop : 310 ! additional loop : 310 311 ! 1) pressure updated 311 312 ! 2) direct redistribution of pressure over the globe … … 339 340 ENDDO 340 341 ENDIF ! (itsub=1) 341 342 342 343 DO ig=1,klon 343 344 ! forecast of frost temperature ztcondsol … … 353 354 ! Condensation or partial sublimation of N2 ice 354 355 if (ztsrf(ig) .LT. ztcondsol(ig)) then ! condensation 355 ! Include a correction to account for the cooling of air near 356 ! Include a correction to account for the cooling of air near 356 357 ! the surface before condensing: 357 358 zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) & … … 359 360 else ! sublimation 360 361 zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) & 361 /(lw_n2*subtimestep) 362 /(lw_n2*subtimestep) 362 363 end if 363 364 … … 369 370 IF((zpicen2(ig)/subtimestep).LE. & 370 371 -zcondices(ig))THEN 371 zcondices(ig) = -zpicen2(ig)/subtimestep 372 zcondices(ig) = -zpicen2(ig)/subtimestep 372 373 pdtsrfc(ig)=(lw_n2/pcapcal(ig))* & 373 374 (zcondices(ig)) … … 376 377 ! Changing N2 ice amount and pressure 377 378 378 pdicen2(ig) = zcondices(ig) 379 pdicen2(ig) = zcondices(ig) 379 380 pdpsrf(ig) = -pdicen2(ig)*g 380 381 ! pdpsrf(ig) = 0. ! OPTION to check impact N2 sub/cond … … 411 412 if(.not.picen2(ig).ge.0.) THEN 412 413 ! if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then 413 print*, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)*ptimestep 414 print*, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)*ptimestep 414 415 ! pdicen2(ig)= -picen2(ig)/ptimestep 415 416 ! else … … 420 421 421 422 ! *************************************************************** 422 ! Correction to account for redistribution between sigma or hybrid 423 ! Correction to account for redistribution between sigma or hybrid 423 424 ! layers when changing surface pressure (and warming/cooling 424 425 ! of the n2 currently changing phase). … … 431 432 ! """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" 432 433 zmflux(1) = -zcondices(ig) 433 DO l=1, nlayer434 DO l=1,klev 434 435 zmflux(l+1) = zmflux(l) -zcondicea(ig,l) & 435 + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1)) 436 ! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld 436 + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1)) 437 ! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld 437 438 if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0. 438 439 END DO 439 440 440 441 ! Mass of each layer 441 ! ------------------ 442 DO l=1, nlayer442 ! ------------------ 443 DO l=1,klev 443 444 masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g 444 445 END DO … … 447 448 ! Corresponding fluxes for T,U,V,Q 448 449 ! """""""""""""""""""""""""""""""" 449 ! averaging operator for TRANSPORT 450 ! averaging operator for TRANSPORT 450 451 ! """""""""""""""""""""""""""""""" 451 452 452 453 ! Subtimestep loop to perform the redistribution gently and simultaneously with 453 ! the other tendencies 454 ! the other tendencies 454 455 ! Estimation of subtimestep (using only the first layer, the most critical) 455 456 dtmax=ptimestep 456 if (zmflux(1).gt.1.e-20) then 457 if (zmflux(1).gt.1.e-20) then 457 458 dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1))) 458 459 endif 459 nsubtimestep= max(nint(ptimestep/dtmax),nint(2.)) 460 nsubtimestep= max(nint(ptimestep/dtmax),nint(2.)) 460 461 subtimestep=ptimestep/float(nsubtimestep) 461 462 462 463 ! New flux for each subtimestep 463 do l=1, nlayer+1464 do l=1,klev+1 464 465 w(l)=-zmflux(l)*subtimestep 465 466 enddo 466 467 ! initializing variables that will vary during subtimestep: 467 do l=1, nlayer468 ztc(l) =pt(ig,l) 469 zu(l) =pu(ig,l) 470 zv(l) =pv(ig,l) 471 do iq=1,nq mx472 zq(l,iq) = pq(ig,l,iq) 468 do l=1,klev 469 ztc(l) =pt(ig,l) 470 zu(l) =pu(ig,l) 471 zv(l) =pv(ig,l) 472 do iq=1,nq 473 zq(l,iq) = pq(ig,l,iq) 473 474 enddo 474 475 end do … … 477 478 ! """""""""""""""""""""" 478 479 do itsub=1,nsubtimestep 479 ! Progressively adding tendancies from other processes. 480 do l=1, nlayer480 ! Progressively adding tendancies from other processes. 481 do l=1,klev 481 482 ztc(l) =ztc(l) +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep 482 483 zu(l) =zu(l) +pdu( ig,l) * subtimestep 483 484 zv(l) =zv(l) +pdv( ig,l) * subtimestep 484 do iq=1,nq mx485 do iq=1,nq 485 486 zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep 486 487 enddo … … 488 489 489 490 ! Change of mass in each layer 490 do l=1, nlayer491 do l=1,klev 491 492 masse(l)=masse(l)+pdpsrf(ig)*subtimestep*(pplev(ig,l) - pplev(ig,l+1))& 492 493 /(g*pplev(ig,1)) 493 494 end do 494 495 495 ztm(2: nlayermx+1)=0.496 zum(2: nlayermx+1)=0.497 zvm(2: nlayermx+1)=0.498 zqm1(1: nlayermx+1)=0.496 ztm(2:klev+1)=0. 497 zum(2:klev+1)=0. 498 zvm(2:klev+1)=0. 499 zqm1(1:klev+1)=0. 499 500 500 501 ! Van Leer scheme: 501 call vl1d( nlayer,ztc,2.,masse,w,ztm)502 call vl1d( nlayer,zu ,2.,masse,w,zum)503 call vl1d( nlayer,zv ,2.,masse,w,zvm)504 do iq=1,nq mx505 do l=1, nlayer502 call vl1d(klev,ztc,2.,masse,w,ztm) 503 call vl1d(klev,zu ,2.,masse,w,zum) 504 call vl1d(klev,zv ,2.,masse,w,zvm) 505 do iq=1,nq 506 do l=1,klev 506 507 zq1(l)=zq(l,iq) 507 508 enddo 508 509 zqm1(1)=zqm(1,iq) 509 call vl1d( nlayer,zq1,2.,masse,w,zqm1)510 do l=2, nlayer510 call vl1d(klev,zq1,2.,masse,w,zqm1) 511 do l=2,klev 511 512 zqm(l,iq)=zqm1(l) 512 513 enddo … … 516 517 if (igcm_n2.ne.0) then 517 518 zqm(1,igcm_n2)=1. 518 do l=1, nlayer-1519 if (w(l)*zqm(l,igcm_n2).gt.zq(l,igcm_n2)*masse(l)) then 519 do l=1,klev-1 520 if (w(l)*zqm(l,igcm_n2).gt.zq(l,igcm_n2)*masse(l)) then 520 521 zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2), & 521 522 (zqm(l,igcm_n2)*w(l) -zq(l,igcm_n2)*masse(l))/w(l+1) ) 522 else 523 else 523 524 exit 524 525 endif … … 528 529 ! Value transfert at the surface interface when condensation sublimation: 529 530 530 if (zmflux(1).lt.0) then 531 if (zmflux(1).lt.0) then 531 532 ! Surface condensation 532 533 zum(1)= zu(1) 533 zvm(1)= zv(1) 534 zvm(1)= zv(1) 534 535 ztm(1) = ztc(1) 535 else 536 else 536 537 ! Surface sublimation: 537 538 ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep 538 zum(1) = 0 539 zvm(1) = 0 539 zum(1) = 0 540 zvm(1) = 0 540 541 end if 541 do iq=1,nq mx542 do iq=1,nq 542 543 zqm(1,iq)=0. ! most tracer do not condense ! 543 544 enddo … … 550 551 DO iq=1,nq 551 552 tname=noms(iq) 552 if (tname(1:4).eq."haze") then 553 if (tname(1:4).eq."haze") then 553 554 !zqm(1,iq)=0.02 554 555 !zqm(1,iq)=-pdicen2(ig)*0.02 … … 559 560 endif 560 561 ENDDO 561 ENDIF 562 ENDIF 563 ztm( nlayer+1)= ztc(nlayer) ! should not be used, but...564 zum( nlayer+1)= zu(nlayer) ! should not be used, but...565 zvm( nlayer+1)= zv(nlayer) ! should not be used, but...566 do iq=1,nq mx567 zqm( nlayer+1,iq)= zq(nlayer,iq)562 ENDIF 563 ENDIF 564 ztm(klev+1)= ztc(klev) ! should not be used, but... 565 zum(klev+1)= zu(klev) ! should not be used, but... 566 zvm(klev+1)= zv(klev) ! should not be used, but... 567 do iq=1,nq 568 zqm(klev+1,iq)= zq(klev,iq) 568 569 enddo 569 570 570 ! Tendencies on T, U, V, Q 571 ! Tendencies on T, U, V, Q 571 572 ! """"""""""""""""""""""" 572 DO l=1, nlayer573 DO l=1,klev 573 574 574 575 ! Tendencies on T … … 591 592 592 593 ! Tendencies on Q 593 do iq=1,nq mx594 do iq=1,nq 594 595 if (iq.eq.igcm_n2) then 595 596 ! SPECIAL Case when the tracer IS N2 : 596 DO l=1, nlayer597 DO l=1,klev 597 598 pdqc(ig,l,iq)= (1/masse(l)) * & 598 599 ( zmflux(l)*(zqm(l,iq) - zq(l,iq)) & … … 601 602 END DO 602 603 else 603 DO l=1, nlayer604 DO l=1,klev 604 605 pdqc(ig,l,iq)= (1/masse(l)) * & 605 606 ( zmflux(l)*(zqm(l,iq) - zq(l,iq)) & … … 610 611 enddo 611 612 ! Update variables at the end of each subtimestep. 612 do l=1, nlayer613 do l=1,klev 613 614 ztc(l) =ztc(l) + zdtsig(ig,l) *subtimestep 614 615 zu(l) =zu(l) + pduc(ig,l) *subtimestep 615 616 zv(l) =zv(l) + pdvc(ig,l) *subtimestep 616 do iq=1,nq mx617 do iq=1,nq 617 618 zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep 618 619 enddo … … 620 621 enddo ! loop on nsubtimestep 621 622 ! Recomputing Total tendencies 622 do l=1, nlayer623 do l=1,klev 623 624 pdtc(ig,l) = (ztc(l) - zt(ig,l) )/ptimestep 624 625 pduc(ig,l) = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep 625 626 pdvc(ig,l) = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep 626 do iq=1,nq mx627 do iq=1,nq 627 628 pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep 628 629 … … 632 633 if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) & 633 634 .lt.0.01) then ! if n2 < 1 % ! 634 pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq) 635 pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq) 635 636 end if 636 637 end if … … 639 640 end do 640 641 ! *******************************TEMPORAIRE ****************** 641 if (klon.eq.1) then 642 if (klon.eq.1) then 642 643 write(*,*) 'nsubtimestep=' ,nsubtimestep 643 644 write(*,*) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g 644 write(*,*) 'masse apres' , masse(1) 645 write(*,*) 'masse apres' , masse(1) 645 646 write(*,*) 'zmflux*DT, l=1' , zmflux(1)*ptimestep 646 647 write(*,*) 'zmflux*DT, l=2' , zmflux(2)*ptimestep … … 653 654 ! *********************************************************** 654 655 end if ! if (condsub) 655 END DO ! loop on ig 656 END DO ! loop on ig 656 657 endif ! not fast 657 658 … … 660 661 if (olkin) then 661 662 DO ig=1,klon 662 if (lati (ig).lt.0) then663 if (latitude(ig).lt.0) then 663 664 pdtsrfc(ig) = max(0.,pdtsrfc(ig)) 664 665 pdpsrf(ig) = 0. 665 666 pdicen2(ig) = 0. 666 do l=1, nlayer667 do l=1,klev 667 668 pdtc(ig,l) = max(0.,zdtlatent(ig,l)) 668 669 pduc(ig,l) = 0. 669 670 pdvc(ig,l) = 0. 670 do iq=1,nq mx671 do iq=1,nq 671 672 pdqc(ig,l,iq) = 0 672 673 enddo … … 681 682 ! *************************************************************** 682 683 683 ! DO l=1, nlayer684 ! DO l=1,klev 684 685 ! DO ig=1,klon 685 686 ! Taux de cond en kg.m-2.pa-1.s-1 … … 689 690 ! END DO 690 691 ! END DO 691 ! call WRITEDIAGFI( ngridmx,'tconda1', &692 ! call WRITEDIAGFI(klon,'tconda1', & 692 693 ! 'Taux de condensation N2 atmospherique /Pa', & 693 694 ! 'kg.m-2.Pa-1.s-1',3,tconda1) 694 ! call WRITEDIAGFI( ngridmx,'tconda2', &695 ! call WRITEDIAGFI(klon,'tconda2', & 695 696 ! 'Taux de condensation N2 atmospherique /m', & 696 697 ! 'kg.m-3.s-1',3,tconda2) … … 698 699 699 700 return 700 end subroutine condens _n2701 end subroutine condense_n2 701 702 702 703 !------------------------------------------------------------------------- … … 729 730 ! tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr)) 730 731 IF (t.ge.35.6) then 731 ! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB 732 ! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB 732 733 ! pcond_n2 = 0.125570*1.e5/vmr*exp((2.5e5*0.98)/296.925*(1./63.1470-1./t)) 733 734 pcond_n2 = 0.125570e5/vmr*exp(825.1241896*(1./63.147-1./t)) … … 743 744 744 745 real function glob_average2d(var) 745 ! Calculates the global average of variable var 746 ! Calculates the global average of variable var 746 747 use comgeomfi_h 748 use dimphy, only: klon 749 USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat 750 use geometry_mod, only: cell_area, latitude 751 747 752 implicit none 748 753 749 754 ! INTEGER klon 750 REAL var( ngridmx)755 REAL var(klon) 751 756 INTEGER ig 752 757 753 758 glob_average2d = 0. 754 DO ig=2, ngridmx-1755 glob_average2d = glob_average2d + var(ig)* area(ig)759 DO ig=2,klon-1 760 glob_average2d = glob_average2d + var(ig)*cell_area(ig) 756 761 END DO 757 glob_average2d = glob_average2d + var(1)* area(1)*iim758 glob_average2d = glob_average2d + var( ngridmx)*area(ngridmx)*iim759 glob_average2d = glob_average2d/(totarea+( area(1)+area(ngridmx))*(iim-1))762 glob_average2d = glob_average2d + var(1)*cell_area(1)*nbp_lon 763 glob_average2d = glob_average2d + var(klon)*cell_area(klon)*nbp_lon 764 glob_average2d = glob_average2d/(totarea+(cell_area(1)+cell_area(klon))*(nbp_lon-1)) 760 765 761 766 end function glob_average2d … … 763 768 ! ***************************************************************** 764 769 765 subroutine vl1d( nlayer,q,pente_max,masse,w,qm)766 ! 770 subroutine vl1d(klev,q,pente_max,masse,w,qm) 771 ! 767 772 ! Operateur de moyenne inter-couche pour calcul de transport type 768 773 ! Van-Leer " pseudo amont " dans la verticale … … 776 781 ! Arguments: 777 782 ! ---------- 778 integer nlayer779 real masse( nlayer),pente_max780 REAL q( nlayer),qm(nlayer+1)781 REAL w( nlayer+1)782 ! 783 ! Local 783 integer klev 784 real masse(klev),pente_max 785 REAL q(klev),qm(klev+1) 786 REAL w(klev+1) 787 ! 788 ! Local 784 789 ! --------- 785 790 ! 786 791 INTEGER l 787 792 ! 788 real dzq( nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax793 real dzq(klev),dzqw(klev),adzqw(klev),dzqmax 789 794 real sigw, Mtot, MQtot 790 integer m 791 792 793 ! On oriente tout dans le sens de la pression 795 integer m 796 797 798 ! On oriente tout dans le sens de la pression 794 799 ! W > 0 WHEN DOWN !!!!!!!!!!!!! 795 800 796 do l=2, nlayer801 do l=2,klev 797 802 dzqw(l)=q(l-1)-q(l) 798 803 adzqw(l)=abs(dzqw(l)) 799 804 enddo 800 805 801 do l=2, nlayer-1806 do l=2,klev-1 802 807 if(dzqw(l)*dzqw(l+1).gt.0.) then 803 808 dzq(l)=0.5*(dzqw(l)+dzqw(l+1)) … … 810 815 811 816 dzq(1)=0. 812 dzq( nlayer)=0.813 814 do l = 1, nlayer-1817 dzq(klev)=0. 818 819 do l = 1,klev-1 815 820 816 821 ! Regular scheme (transfered mass < layer mass) … … 829 834 Mtot = masse(m) 830 835 MQtot = masse(m)*q(m) 831 !if (m.lt. nlayer) then ! because some compilers will have problems832 ! ! evaluating masse( nlayer+1)833 do while ((m.lt. nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))836 !if (m.lt.klev) then ! because some compilers will have problems 837 ! ! evaluating masse(klev+1) 838 do while ((m.lt.klev).and.(w(l+1).gt.(Mtot+masse(m+1)))) 834 839 m=m+1 835 840 Mtot = Mtot + masse(m) 836 841 MQtot = MQtot + masse(m)*q(m) 837 ! if (m.eq. nlayer) exit842 ! if (m.eq.klev) exit 838 843 end do 839 844 !endif 840 if (m.lt. nlayer) then845 if (m.lt.klev) then 841 846 sigw=(w(l+1)-Mtot)/masse(m+1) 842 847 qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)* & … … 848 853 stop 849 854 end if 850 else ! if(w(l+1).lt.0) 851 m = l-1 855 else ! if(w(l+1).lt.0) 856 m = l-1 852 857 Mtot = masse(m+1) 853 858 MQtot = masse(m+1)*q(m+1) … … 863 868 if (m.gt.0) then 864 869 sigw=(w(l+1)+Mtot)/masse(m) 865 qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* & 870 qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* & 866 871 (q(m)-0.5*(1.+sigw)*dzq(m)) ) 867 872 else 868 873 qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1)) 869 end if 874 end if 870 875 end if 871 876 enddo -
trunk/LMDZ.PLUTO/libf/phypluto/convadj.F
r3184 r3195 7 7 USE tracer_h 8 8 use comcstfi_mod, only: g 9 use callkeys_mod, only: tracer ,water9 use callkeys_mod, only: tracer 10 10 11 11 implicit none 12 12 13 13 !================================================================== 14 ! 14 ! 15 15 ! Purpose 16 16 ! ------- … … 23 23 ! Original author unknown. 24 24 ! Modif. 2005 by F. Forget. 25 ! 25 ! 26 26 !================================================================== 27 27 … … 84 84 ENDDO 85 85 86 if(tracer) then 86 if(tracer) then 87 87 DO iq =1, nq 88 88 DO l=1,nlay … … 118 118 ENDDO 119 119 ENDDO 120 120 121 121 ! Make a list of them 122 122 jcnt=0 … … 136 136 137 137 i = jadrs(jj) 138 138 139 139 ! Calculate sigma in this column 140 140 DO l=1,nlay+1 141 141 sig(l)=pplev(i,l)/pplev(i,1) 142 142 143 143 ENDDO 144 144 DO l=1,nlay … … 154 154 IF (l2 .GT. nlay) EXIT 155 155 IF (zhc(i, l2).LT.zhc(i, l2-1)) THEN 156 156 157 157 ! l2 is the highest level of the unstable column 158 158 159 159 l1 = l2 - 1 160 160 l = l1 … … 170 170 zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm 171 171 zhmc = zhm 172 172 173 173 ! do we have to extend the column downwards? 174 174 175 175 down = .false. 176 176 IF (l1 .ne. 1) then !-- and then … … 179 179 END IF 180 180 END IF 181 181 182 182 ! this could be a problem... 183 183 184 184 if (down) then 185 185 186 186 l1 = l1 - 1 187 187 l = l1 188 188 189 189 else 190 190 191 191 ! can we extend the column upwards? 192 192 193 193 if (l2 .eq. nlay) exit 194 194 195 195 if (zhc(i, l2+1) .ge. zhmc) exit 196 196 … … 224 224 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 225 225 226 ! to conserve tracers/ KE, we must calculate zum, zvm and zqm using 227 ! the up-to-date column values. If we do not do this, there are cases 226 ! to conserve tracers/ KE, we must calculate zum, zvm and zqm using 227 ! the up-to-date column values. If we do not do this, there are cases 228 228 ! where convection stops at one level and starts at the next where we 229 229 ! can break conservation of stuff (particularly tracers) significantly. … … 249 249 zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l)) 250 250 do iq=1,nq 251 ! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq)) 251 ! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq)) 252 252 zq2(i,l,iq)=zqm(iq) 253 253 end do … … 266 266 ! check conservation 267 267 cadjncons=0.0 268 if(water)then269 do l = 1, nlay270 masse = (pplev(i,l) - pplev(i,l+1))/g271 iq = igcm_h2o_vap272 cadjncons = cadjncons +273 & masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep274 end do275 endif276 268 277 269 if(cadjncons.lt.-1.e-6)then … … 283 275 do l = 1, nlay 284 276 print*,'dsig = ',dsig(l) 285 end do 277 end do 286 278 do l = 1, nlay 287 279 print*,'dsig = ',dsig(l) … … 311 303 print*,'zh2(ig,:) = ',zh2(i,l) 312 304 end do 313 do l = 1, nlay314 print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_vap)315 end do316 do l = 1, nlay317 print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_vap)318 end do319 print*,'zqm(vap) = ',zqm(igcm_h2o_vap)305 ! do l = 1, nlay 306 ! print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_vap) 307 ! end do 308 ! do l = 1, nlay 309 ! print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_vap) 310 ! end do 311 ! print*,'zqm(vap) = ',zqm(igcm_h2o_vap) 320 312 print*,'jadrs=',jadrs 321 313 … … 335 327 ENDDO 336 328 337 if(tracer) then 329 if(tracer) then 338 330 do iq=1, nq 339 331 do l=1,nlay 340 332 DO ig=1,ngrid 341 pdqadj(ig,l,iq)=(zq2(ig,l,iq)-zq(ig,l,iq))/ptimestep 342 end do 343 end do 344 end do 333 pdqadj(ig,l,iq)=(zq2(ig,l,iq)-zq(ig,l,iq))/ptimestep 334 end do 335 end do 336 end do 345 337 end if 346 338 … … 349 341 ! if (ngrid.eq.1) then 350 342 ! ig=1 351 ! iq =1 352 ! write(*,*)'**** l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)' 343 ! iq =1 344 ! write(*,*)'**** l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)' 353 345 ! do l=nlay,1,-1 354 346 ! write(*,*) l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq) -
trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90
r3193 r3195 24 24 use planetwide_mod, only: planetwide_sumval 25 25 use callkeys_mod 26 use surfdat_h 26 27 use wstats_mod, only: callstats 27 28 use ioipsl_getin_p_mod, only : getin_p … … 277 278 if (is_master) write(*,*) trim(rname)//": calladj = ",calladj 278 279 279 if (is_master) write(*,*) trim(rname)//": call N2 condensation ?"280 n2cond=.false. ! default value281 call getin_p("n2cond",n2cond)282 if (is_master) write(*,*) trim(rname)//": n2cond = ",n2cond283 ! Test of incompatibility284 if (n2cond.and.(.not.tracer)) then285 call abort_physic(rname,"Error: We need a N2 ice tracer to condense N2!",1)286 endif287 288 if (is_master) write(*,*) trim(rname)//": N2 supersaturation level ?"289 n2supsat=1.0 ! default value290 call getin_p("n2supsat",n2supsat)291 if (is_master) write(*,*) trim(rname)//": n2supsat = ",n2supsat292 293 280 if (is_master) write(*,*) trim(rname)//& 294 281 ": Radiative timescale for Newtonian cooling (in Earth days)?" … … 479 466 endif 480 467 481 if (is_master) write(*,*)trim(rname)//& 482 ": Number mixing ratio of N2 ice particles:" 483 Nmix_n2=1.e6 ! default value 484 call getin_p("Nmix_n2",Nmix_n2) 485 if (is_master) write(*,*)trim(rname)//": Nmix_n2 = ",Nmix_n2 486 487 if (is_master) write(*,*)trim(rname)//& 468 if (is_master) write(*,*)trim(rname)//& 488 469 "Number of radiatively active aerosols:" 489 470 naerkind=0 ! default value … … 491 472 if (is_master) write(*,*)trim(rname)//": naerkind = ",naerkind 492 473 493 if (is_master) write(*,*)trim(rname)//": Opacity of dust (if used):" 494 dusttau=0. ! default value 495 call getin_p("dusttau",dusttau) 496 if (is_master) write(*,*)trim(rname)//": dusttau = ",dusttau 497 498 if (is_master) write(*,*)trim(rname)//": Radiatively active N2 aerosols?" 499 aeron2=.false. ! default value 500 call getin_p("aeron2",aeron2) 501 if (is_master) write(*,*)trim(rname)//": aeron2 = ",aeron2 502 503 if (is_master) write(*,*)trim(rname)//": Fixed N2 aerosol distribution?" 504 aerofixn2=.false. ! default value 505 call getin_p("aerofixn2",aerofixn2) 506 if (is_master) write(*,*)trim(rname)//": aerofixn2 = ",aerofixn2 507 508 if (is_master) write(*,*)trim(rname)//& 509 ": Radiatively active sulfuric acid aerosols?" 510 aeroh2so4=.false. ! default value 511 call getin_p("aeroh2so4",aeroh2so4) 512 if (is_master) write(*,*)trim(rname)//": aeroh2so4 = ",aeroh2so4 513 514 if (is_master) write(*,*)trim(rname)//& 515 ": Radiatively active auroral aerosols?" 516 aeroaurora=.false. ! default value 517 call getin_p("aeroaurora",aeroaurora) 518 if (is_master) write(*,*)trim(rname)//": aeroaurora = ",aeroaurora 519 520 if (is_master) write(*,*)trim(rname)//& 521 ": Radiatively active Generic Condensable aerosols?" 522 aerogeneric=0 ! default value 523 call getin_p("aerogeneric",aerogeneric) 524 if (is_master) write(*,*)trim(rname)//":aerogeneric",aerogeneric 525 474 475 !*************************************************************** 476 !! TRACERS options 477 478 if (is_master)write(*,*)trim(rname)//& 479 "call N2 condensation ?" 480 N2cond=.true. ! default value 481 call getin_p("N2cond",N2cond) 482 if (is_master)write(*,*)trim(rname)//& 483 " N2cond = ",N2cond 484 485 if (is_master)write(*,*)trim(rname)//& 486 "call N2 cloud condensation ?" 487 condensn2=.false. ! default value 488 call getin_p("condensn2",condensn2) 489 if (is_master)write(*,*)trim(rname)//& 490 "condensn2 = ",condensn2 491 492 if (is_master)write(*,*)trim(rname)//& 493 "call no N2 frost formation?" 494 no_n2frost=.false. ! default value 495 call getin_p("no_n2frost",no_n2frost) 496 if (is_master)write(*,*)trim(rname)//& 497 "no_n2frost = ",no_n2frost 498 499 if (is_master)write(*,*)trim(rname)//& 500 "N2 condensation subtimestep?" 501 nbsub=20 ! default value 502 call getin_p("nbsub",nbsub) 503 if (is_master)write(*,*)trim(rname)//& 504 " nbsub = ",nbsub 505 506 if (is_master)write(*,*)trim(rname)//& 507 "Gravitationnal sedimentation ?" 508 sedimentation=.true. ! default value 509 call getin_p("sedimentation",sedimentation) 510 if (is_master)write(*,*)trim(rname)//& 511 " sedimentation = ",sedimentation 512 513 if (is_master)write(*,*)trim(rname)//& 514 "Compute glacial flow ?" 515 glaflow=.false. ! default value 516 call getin_p("glaflow",glaflow) 517 if (is_master)write(*,*)trim(rname)//& 518 " glaflow = ",glaflow 519 520 if (is_master)write(*,*)trim(rname)//& 521 "Compute methane cycle ?" 522 methane=.false. ! default value 523 call getin_p("methane",methane) 524 if (is_master)write(*,*)trim(rname)//& 525 " methane = ",methane 526 condmetsurf=.true. ! default value 527 call getin_p("condmetsurf",condmetsurf) 528 if (is_master)write(*,*)trim(rname)//& 529 " condmetsurf = ",condmetsurf 530 531 if (is_master)write(*,*)trim(rname)//& 532 "Compute methane clouds ?" 533 metcloud=.false. ! default value 534 call getin_p("metcloud",metcloud) 535 if (is_master)write(*,*)trim(rname)//& 536 " metcloud = ",metcloud 537 538 if (is_master)write(*,*)trim(rname)//& 539 "Compute CO cycle ?" 540 carbox=.false. ! default value 541 call getin_p("carbox",carbox) 542 if (is_master)write(*,*)trim(rname)//& 543 " carbox = ",carbox 544 condcosurf=.true. ! default value 545 call getin_p("condcosurf",condcosurf) 546 if (is_master)write(*,*)trim(rname)//& 547 " condcosurf = ",condcosurf 548 549 if (is_master)write(*,*)trim(rname)//& 550 "Compute CO clouds ?" 551 monoxcloud=.false. ! default value 552 call getin_p("monoxcloud",monoxcloud) 553 if (is_master)write(*,*)trim(rname)//& 554 " monoxcloud = ",monoxcloud 555 556 if (is_master)write(*,*)trim(rname)//& 557 "atmospheric redistribution (s):" 558 tau_n2=1. ! default value 559 call getin_p("tau_n2",tau_n2) 560 if (is_master)write(*,*)trim(rname)//& 561 " tau_n2 = ",tau_n2 562 tau_ch4=1.E7 ! default value 563 call getin_p("tau_ch4",tau_ch4) 564 if (is_master)write(*,*)trim(rname)//& 565 " tau_ch4 = ",tau_ch4 566 tau_co=1. ! default value 567 call getin_p("tau_co",tau_co) 568 if (is_master)write(*,*)trim(rname)//& 569 " tau_co = ",tau_co 570 tau_prechaze=1. ! default value 571 call getin_p("tau_prechaze",tau_prechaze) 572 if (is_master)write(*,*)trim(rname)//& 573 " tau_prechaze = ",tau_prechaze 574 575 if (is_master)write(*,*)trim(rname)//& 576 "Day fraction for limited cold trap in SP?" 577 dayfrac=0. ! default value 578 call getin_p("dayfrac",dayfrac) 579 if (is_master)write(*,*)trim(rname)//& 580 " dayfrac = ",dayfrac 581 thresh_non2=0. ! default value 582 call getin_p("thresh_non2",thresh_non2) 583 if (is_master)write(*,*)trim(rname)//& 584 " thresh_non2 = ",thresh_non2 585 586 !*************************************************************** 587 !! Haze options 588 589 if (is_master)write(*,*)trim(rname)//& 590 "Production of haze ?" 591 haze=.false. ! default value 592 call getin_p("haze",haze) 593 if (is_master)write(*,*)trim(rname)//& 594 " haze = ",haze 595 hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed 596 call getin_p("hazeconservch4",hazeconservch4) 597 if (is_master)write(*,*)trim(rname)//& 598 "hazconservch4 = ",hazeconservch4 599 if (is_master)write(*,*)trim(rname)//& 600 "Production of haze (fast version) ?" 601 fasthaze=.false. ! default value 602 call getin_p("fasthaze",fasthaze) 603 if (is_master)write(*,*)trim(rname)//& 604 "fasthaze = ",fasthaze 605 606 if (is_master)write(*,*)trim(rname)//& 607 "Add source of haze ?" 608 source_haze=.false. ! default value 609 call getin_p("source_haze",source_haze) 610 if (is_master)write(*,*)trim(rname)//& 611 " source_haze = ",source_haze 612 mode_hs=0 ! mode haze source default value 613 call getin_p("mode_hs",mode_hs) 614 if (is_master)write(*,*)trim(rname)//& 615 " mode_hs",mode_hs 616 kfix=1 ! default value 617 call getin_p("kfix",kfix) 618 if (is_master)write(*,*)trim(rname)//& 619 "levels kfix",kfix 620 fracsource=0.1 ! default value 621 call getin_p("fracsource",fracsource) 622 if (is_master)write(*,*)trim(rname)//& 623 " fracsource",fracsource 624 latsource=30. ! default value 625 call getin_p("latsource",latsource) 626 if (is_master)write(*,*)trim(rname)//& 627 " latsource",latsource 628 lonsource=180. ! default value 629 call getin_p("lonsource",lonsource) 630 if (is_master)write(*,*)trim(rname)//& 631 " lonsource",lonsource 632 633 if (is_master)write(*,*)trim(rname)//& 634 "Radiatively active haze ?" 635 aerohaze=.false. ! default value 636 call getin_p("aerohaze",aerohaze) 637 if (is_master)write(*,*)trim(rname)//& 638 "aerohaze = ",aerohaze 639 640 if (is_master)write(*,*)trim(rname)//& 641 "Haze monomer radius ?" 642 rad_haze=10.e-9 ! default value 643 call getin_p("rad_haze",rad_haze) 644 if (is_master)write(*,*)trim(rname)//& 645 "rad_haze = ",rad_haze 646 647 if (is_master)write(*,*)trim(rname)//& 648 "fractal particle ?" 649 fractal=.false. ! default value 650 call getin_p("fractal",fractal) 651 if (is_master)write(*,*)trim(rname)//& 652 "fractal = ",fractal 653 nb_monomer=10 ! default value 654 call getin_p("nb_monomer",nb_monomer) 655 if (is_master)write(*,*)trim(rname)//& 656 "nb_monomer = ",nb_monomer 657 658 if (is_master)write(*,*)trim(rname)//& 659 "fixed radius profile from txt file ?" 660 haze_radproffix=.false. ! default value 661 call getin_p("haze_radproffix",haze_radproffix) 662 if (is_master)write(*,*)trim(rname)//& 663 "haze_radproffix = ",haze_radproffix 664 if (is_master)write(*,*)trim(rname)//& 665 "fixed MMR profile from txt file ?" 666 haze_proffix=.false. ! default value 667 call getin_p("haze_proffix",haze_proffix) 668 if (is_master)write(*,*)trim(rname)//& 669 "haze_proffix = ",haze_proffix 670 671 if (is_master)write(*,*)trim(rname)//& 672 "Number mix ratio of haze particles for co clouds:" 673 Nmix_co=100000. ! default value 674 call getin_p("Nmix_co",Nmix_co) 675 if (is_master)write(*,*)trim(rname)//& 676 " Nmix_co = ",Nmix_co 677 678 if (is_master)write(*,*)trim(rname)//& 679 "Number mix ratio of haze particles for ch4 clouds:" 680 Nmix_ch4=100000. ! default value 681 call getin_p("Nmix_ch4",Nmix_ch4) 682 if (is_master)write(*,*)trim(rname)//& 683 " Nmix_ch4 = ",Nmix_ch4 684 685 !*************************************************************** 686 !! Surface properties 687 688 !*********** N2 ********************************* 689 690 if (is_master)write(*,*)trim(rname)//& 691 "Mode for changing N2 albedo" 692 mode_n2=0 ! default value 693 call getin_p("mode_n2",mode_n2) 694 if (is_master)write(*,*)trim(rname)//& 695 " mode_n2 = ",mode_n2 696 thres_n2ice=1. ! default value 697 call getin_p("thres_n2ice",thres_n2ice) 698 if (is_master)write(*,*)trim(rname)//& 699 " thres_n2ice = ",thres_n2ice 700 701 if (is_master)write(*,*)trim(rname)//& 702 "Diff of N2 albedo with thickness" 703 deltab=0. ! default value 704 call getin_p("deltab",deltab) 705 if (is_master)write(*,*)trim(rname)//& 706 " deltab = ",deltab 707 708 if (is_master)write(*,*)trim(rname)//& 709 "albedo N2 beta " 710 alb_n2b=0.7 ! default value 711 call getin_p("alb_n2b",alb_n2b) 712 if (is_master)write(*,*)trim(rname)//& 713 " alb_n2b = ",alb_n2b 714 715 if (is_master)write(*,*)trim(rname)//& 716 "albedo N2 alpha " 717 alb_n2a=0.7 ! default value 718 call getin_p("alb_n2a",alb_n2a) 719 if (is_master)write(*,*)trim(rname)//& 720 " alb_n2a = ",alb_n2a 721 722 if (is_master)write(*,*)trim(rname)//& 723 "emis N2 beta " 724 emis_n2b=0.7 ! default value 725 call getin_p("emis_n2b",emis_n2b) 726 if (is_master)write(*,*)trim(rname)//& 727 " emis_n2b = ",emis_n2b 728 729 if (is_master)write(*,*)trim(rname)//& 730 "emis N2 alpha " 731 emis_n2a=0.7 ! default value 732 call getin_p("emis_n2a",emis_n2a) 733 if (is_master)write(*,*)trim(rname)//& 734 " emis_n2a = ",emis_n2a 735 736 !*********** CH4 ********************************* 737 738 if (is_master)write(*,*)trim(rname)//& 739 "Mode for changing CH4 albedo" 740 mode_ch4=0 ! default value 741 call getin_p("mode_ch4",mode_ch4) 742 if (is_master)write(*,*)trim(rname)//& 743 " mode_ch4 = ",mode_ch4 744 feedback_met=0 ! default value 745 call getin_p("feedback_met",feedback_met) 746 if (is_master)write(*,*)trim(rname)//& 747 " feedback_met = ",feedback_met 748 thres_ch4ice=1. ! default value, mm 749 call getin_p("thres_ch4ice",thres_ch4ice) 750 if (is_master)write(*,*)trim(rname)//& 751 " thres_ch4ice = ",thres_ch4ice 752 fdch4_latn=-1. 753 fdch4_lats=0. 754 fdch4_lone=0. 755 fdch4_lonw=-1. 756 fdch4_depalb=0.5 757 fdch4_finalb=0.95 758 fdch4_maxalb=0.99 759 fdch4_ampl=200. 760 fdch4_maxice=100. 761 call getin_p(" fdch4_latn",fdch4_latn) 762 call getin_p(" fdch4_lats",fdch4_lats) 763 call getin_p(" fdch4_lone",fdch4_lone) 764 call getin_p(" fdch4_lonw",fdch4_lonw) 765 call getin_p(" fdch4_depalb",fdch4_depalb) 766 call getin_p(" fdch4_finalb",fdch4_finalb) 767 call getin_p(" fdch4_maxalb",fdch4_maxalb) 768 call getin_p(" fdch4_maxice",fdch4_maxice) 769 call getin_p(" fdch4_ampl",fdch4_ampl) 770 if (is_master)write(*,*)trim(rname)//& 771 " Values for albedo feedback = ",fdch4_latn,& 772 fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,& 773 fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl 774 775 if (is_master)write(*,*)trim(rname)//& 776 "Latitude for diff albedo methane" 777 metlateq=25. ! default value 778 call getin_p("metlateq",metlateq) 779 if (is_master)write(*,*)trim(rname)//& 780 " metlateq = ",metlateq 781 782 if (is_master)write(*,*)trim(rname)//& 783 "Ls1 and Ls2 for change of ch4 albedo" 784 metls1=-1. ! default value 785 metls2=-2. ! default value 786 call getin_p("metls1",metls1) 787 call getin_p("metls2",metls2) 788 789 if (is_master)write(*,*)trim(rname)//& 790 "albedo CH4 " 791 alb_ch4=0.5 ! default value 792 call getin_p("alb_ch4",alb_ch4) 793 if (is_master)write(*,*)trim(rname)//& 794 " alb_ch4 = ",alb_ch4 795 796 if (is_master)write(*,*)trim(rname)//& 797 "albedo equatorial CH4 " 798 alb_ch4_eq=alb_ch4 ! default value 799 call getin_p("alb_ch4_eq",alb_ch4_eq) 800 if (is_master)write(*,*)trim(rname)//& 801 " alb_ch4_eq = ",alb_ch4_eq 802 803 if (is_master)write(*,*)trim(rname)//& 804 "albedo south hemis CH4 " 805 alb_ch4_s=alb_ch4 ! default value 806 call getin_p("alb_ch4_s",alb_ch4_s) 807 if (is_master)write(*,*)trim(rname)//& 808 " alb_ch4_s = ",alb_ch4_s 809 810 if (is_master)write(*,*)trim(rname)//& 811 "emis CH4 " 812 emis_ch4=0.5 ! default value 813 call getin_p("emis_ch4",emis_ch4) 814 if (is_master)write(*,*)trim(rname)//& 815 " emis_ch4 = ",emis_ch4 816 817 if (is_master)write(*,*)trim(rname)//& 818 "CH4 lag for n2 sublimation limitation" 819 ch4lag=.false. ! default value 820 latlag=-90. ! default value 821 vmrlag=1. ! default value 822 call getin_p("ch4lag",ch4lag) 823 call getin_p("latlag",latlag) 824 call getin_p("vmrlag",vmrlag) 825 if (is_master)write(*,*)trim(rname)//& 826 " ch4lag = ",ch4lag 827 if (is_master)write(*,*)trim(rname)//& 828 " latlag = ",latlag 829 if (is_master)write(*,*)trim(rname)//& 830 " vmrlag = ",vmrlag 831 832 if (is_master)write(*,*)trim(rname)//& 833 "Max temperature for surface ?" 834 tsurfmax=.false. ! default value 835 albmin_ch4=0.3 ! default value 836 call getin_p("tsurfmax",tsurfmax) 837 call getin_p("albmin_ch4",albmin_ch4) 838 if (is_master)write(*,*)trim(rname)//& 839 " tsurfmax = ",tsurfmax 840 if (is_master)write(*,*)trim(rname)//& 841 " albmin_ch4 = ",albmin_ch4 842 843 !*********** CO ********************************* 844 845 if (is_master)write(*,*)trim(rname)//& 846 "albedo CO " 847 alb_co=0.4 ! default value 848 call getin_p("alb_co",alb_co) 849 if (is_master)write(*,*)trim(rname)//& 850 " alb_co = ",alb_co 851 thres_coice=1. ! default value, mm 852 call getin_p("thres_coice",thres_coice) 853 if (is_master)write(*,*)trim(rname)//& 854 " thres_coice = ",thres_coice 855 856 if (is_master)write(*,*)trim(rname)//& 857 "emis CO " 858 emis_co=0.4 ! default value 859 call getin_p("emis_co",emis_co) 860 if (is_master)write(*,*)trim(rname)//& 861 " emis_co = ",emis_co 862 863 !*********** THOLINS ********************************* 864 if (is_master)write(*,*)trim(rname)//& 865 "Mode for changing tholins albedo/emis" 866 mode_tholins=0 ! default value 867 call getin_p("mode_tholins",mode_tholins) 868 if (is_master)write(*,*)trim(rname)//& 869 " mode_tholins = ",mode_tholins 870 871 if (is_master)write(*,*)trim(rname)//& 872 "albedo tho " 873 alb_tho=0.1 ! default value 874 call getin_p("alb_tho",alb_tho) 875 if (is_master)write(*,*)trim(rname)//& 876 " alb_tho = ",alb_tho 877 878 if (is_master)write(*,*)trim(rname)//& 879 "albedo tho eq" 880 alb_tho_eq=0.1 ! default value 881 call getin_p("alb_tho_eq",alb_tho_eq) 882 if (is_master)write(*,*)trim(rname)//& 883 " alb_tho_eq = ",alb_tho_eq 884 885 if (is_master)write(*,*)trim(rname)//& 886 "emis tholins " 887 emis_tho=1. ! default value 888 call getin_p("emis_tho",emis_tho) 889 if (is_master)write(*,*)trim(rname)//& 890 " emis_tho = ",emis_tho 891 892 if (is_master)write(*,*)trim(rname)//& 893 "emis tholins eq" 894 emis_tho_eq=1. ! default value 895 call getin_p("emis_tho_eq",emis_tho_eq) 896 if (is_master)write(*,*)trim(rname)//& 897 " emis_tho_eq = ",emis_tho_eq 898 899 if (is_master)write(*,*)trim(rname)//& 900 "Latitude for diff albedo tholins" 901 tholateq=25. ! default value 902 call getin_p("tholateq",tholateq) 903 if (is_master)write(*,*)trim(rname)//& 904 " tholateq = ",tholateq 905 tholatn=-1. 906 tholats=0. 907 tholone=0. 908 tholonw=-1. 909 alb_tho_spe=0.1 ! default value 910 emis_tho_spe=1. ! default value 911 call getin_p(" tholatn",tholatn) 912 call getin_p(" tholats",tholats) 913 call getin_p(" tholone",tholone) 914 call getin_p(" tholonw",tholonw) 915 if (is_master)write(*,*)trim(rname)//& 916 " Values for special tholins albedo = ",tholatn,& 917 tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe 918 919 if (is_master)write(*,*)trim(rname)//& 920 "Specific albedo" 921 spelon1=-180. ! default value 922 spelon2=180. ! default value 923 spelat1=-90. ! default value 924 spelat2=90. ! default value 925 specalb=.false. ! default value 926 if (is_master)write(*,*)trim(rname)//& 927 "albedo/emis spe " 928 albspe=0.1 ! default value 929 emispe=1. ! default value 930 call getin_p("spelon1",spelon1) 931 call getin_p("spelon2",spelon2) 932 call getin_p("spelat1",spelat1) 933 call getin_p("spelat2",spelat2) 934 call getin_p("specalb",specalb) 935 call getin_p("albspe",albspe) 936 call getin_p("emispe",emispe) 937 938 if (is_master)write(*,*)trim(rname)//& 939 " specific = ",specalb 940 941 !********************** TI ***************************** 942 943 if (is_master)write(*,*)trim(rname)//& 944 "Change TI with time" 945 changeti=.false. ! default value 946 call getin_p("changeti",changeti) 947 if (is_master)write(*,*)trim(rname)//& 948 " changeti = ",changeti 949 changetid=.false. ! default value for diurnal TI 950 call getin_p("changetid",changetid) 951 if (is_master)write(*,*)trim(rname)//& 952 " changetid = ",changetid 953 954 if (is_master)write(*,*)trim(rname)//& 955 "IT N2 " 956 ITN2=800. ! default value 957 call getin_p("ITN2",ITN2) 958 if (is_master)write(*,*)trim(rname)//& 959 " ITN2 = ",ITN2 960 ITN2d=20. ! default value 961 call getin_p("ITN2d",ITN2d) 962 if (is_master)write(*,*)trim(rname)//& 963 " ITN2d = ",ITN2d 964 965 if (is_master)write(*,*)trim(rname)//& 966 "IT CH4" 967 ITCH4=800. ! default value 968 call getin_p("ITCH4",ITCH4) 969 if (is_master)write(*,*)trim(rname)//& 970 " ITCH4 = ",ITCH4 971 ITCH4d=20. ! default value 972 call getin_p("ITCH4d",ITCH4d) 973 if (is_master)write(*,*)trim(rname)//& 974 " ITCH4d = ",ITCH4d 975 976 if (is_master)write(*,*)trim(rname)//& 977 "IT H2O" 978 ITH2O=800. ! default value 979 call getin_p("ITH2O",ITH2O) 980 if (is_master)write(*,*)trim(rname)//& 981 " ITH2O = ",ITH2O 982 ITH2Od=20. ! default value 983 call getin_p("ITH2Od",ITH2Od) 984 if (is_master)write(*,*)trim(rname)//& 985 " ITH2Od = ",ITH2Od 526 986 527 987 !================================= -
trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90
r3193 r3195 3 3 use surfdat_h, ONLY: dryness 4 4 USE tracer_h 5 USE callkeys_mod, only: aerohaze,nb_monomer,haze,fractal,fasthaze,rad_haze 5 6 USE recombin_corrk_mod, ONLY: ini_recombin 6 7 USE mod_phys_lmdz_para, only: is_master, bcast 7 8 use generic_cloud_common_h 9 use aerosol_mod, only: iaero_haze,i_haze 8 10 IMPLICIT NONE 9 11 !======================================================================= 10 12 ! subject: 11 13 ! -------- 12 ! Initialization related to tracer 14 ! Initialization related to tracer 13 15 ! (transported dust, water, chemical species, ice...) 14 16 ! … … 23 25 ! Ehouarn Millour (oct. 2008) identify tracers by their names 24 26 ! Y Jaziri & J. Vatant d'Ollone (2020) : Modern traceur.def 25 ! L Teinturier (2022): Tracer names are now read here instead of 26 ! inside interfaces 27 ! L Teinturier (2022): Tracer names are now read here instead of 28 ! inside interfaces 27 29 !======================================================================= 28 30 … … 30 32 31 33 character(len=500) :: tracline ! to read traceur.def lines 32 ! character(len=30) :: txt ! to store some text33 34 integer :: blank !to store the index of 1st blank when reading tracers names 34 35 integer iq,ig,count,ierr 36 real r0_lift , reff_lift, rho_haze 37 integer nqhaze(nq) ! to store haze tracers 38 integer i, ia, block 39 character(len=20) :: txt ! to store some text 40 CHARACTER(LEN=20) :: tracername ! to temporarily store text 35 41 36 42 !----------------------------------------------------------------------- … … 86 92 ENDDO 87 93 ENDIF ! if modern or standard traceur.def 88 94 89 95 endif ! of if (is_master) 90 96 91 97 ! share the information with other procs/threads (if any) 92 98 CALL bcast(nqtot) 93 99 CALL bcast(moderntracdef) 94 100 95 101 !! ----------------------------------------------------------------------- 96 102 !! For the moment number of tracers in dynamics and physics are equal … … 108 114 IF (.NOT.ALLOCATED(alpha_devil)) ALLOCATE(alpha_devil(nq)) 109 115 IF (.NOT.ALLOCATED(qextrhor)) ALLOCATE(qextrhor(nq)) 110 IF (.NOT.ALLOCATED(igcm_dustbin)) ALLOCATE(igcm_dustbin(nq))116 ! IF (.NOT.ALLOCATED(igcm_dustbin)) ALLOCATE(igcm_dustbin(nq)) 111 117 IF (.NOT.ALLOCATED(is_chim)) ALLOCATE(is_chim(nqtot)) 112 118 IF (.NOT.ALLOCATED(is_rad)) ALLOCATE(is_rad(nqtot)) … … 140 146 is_recomb_qset(:) = 0 141 147 is_recomb_qotf(:) = 0 142 148 143 149 ! Added by JVO 2017 : these arrays are handled later 144 150 ! -> initialization is the least we can do, please !!! … … 163 169 164 170 ! Initialization: Read tracers names from traceur.def 165 do iq=1,nq 171 do iq=1,nq 166 172 if (is_master) read(407,'(A)') tracline 167 173 call bcast(tracline) … … 173 179 ! 0. initialize tracer indexes to zero: 174 180 ! NB: igcm_* indexes are commons in 'tracer.h' 175 do iq=1,nq176 igcm_dustbin(iq)=0177 enddo178 igcm_dust_mass=0179 igcm_dust_number=0180 igcm_h2o_vap=0181 igcm_h2o_ice=0182 igcm_co2=0183 igcm_co=0184 igcm_o=0185 igcm_o1d=0186 igcm_o2=0187 igcm_o3=0188 igcm_h=0189 igcm_h2=0190 igcm_oh=0191 igcm_ho2=0192 igcm_h2o2=0193 181 igcm_n2=0 194 igcm_n=0 195 igcm_n2d=0 196 igcm_no=0 197 igcm_no2=0 198 igcm_ar=0 199 igcm_ar_n2=0 200 igcm_n2_ice=0 201 202 igcm_ch4=0 203 igcm_ch3=0 204 igcm_ch=0 205 igcm_3ch2=0 206 igcm_1ch2=0 207 igcm_cho=0 208 igcm_ch2o=0 209 igcm_ch3o=0 210 igcm_c=0 211 igcm_c2=0 212 igcm_c2h=0 213 igcm_c2h2=0 214 igcm_c2h3=0 215 igcm_c2h4=0 216 igcm_c2h6=0 217 igcm_ch2co=0 218 igcm_ch3co=0 219 igcm_hcaer=0 182 igcm_ch4_gas=0 183 igcm_ch4_ice=0 184 igcm_prec_haze=0 185 igcm_co_gas=0 186 igcm_co_ice=0 187 188 nqhaze(:)=0 189 i=0 190 DO iq=1,nq 191 txt=noms(iq) 192 IF (txt(1:4).eq."haze") THEN 193 i=i+1 194 nqhaze(i)=iq 195 ENDIF 196 ENDDO 197 if ((haze.or.fasthaze).and.i==0) then 198 print*, 'Haze active but no haze tracer in traceur.def' 199 stop 200 endif 201 202 igcm_haze=0 203 igcm_haze10=0 204 igcm_haze30=0 205 igcm_haze50=0 206 igcm_haze100=0 207 208 ! Eddy diffusion tracers 209 igcm_eddy1e6=0 210 igcm_eddy1e7=0 211 igcm_eddy5e7=0 212 igcm_eddy1e8=0 213 igcm_eddy5e8=0 214 write(*,*) 'initracer: noms() ', noms 215 rho_n2=1000 ! n2 ice 216 rho_ch4_ice=520. ! rho ch4 ice 217 rho_co_ice=520. ! rho ch4 ice 218 rho_haze=800. ! haze 220 219 221 220 ! 1. find dust tracers … … 224 223 ! 2. find chemistry and water tracers 225 224 do iq=1,nq 226 if (noms(iq).eq."co2") then227 igcm_co2=iq228 mmol(igcm_co2)=44.229 count=count+1230 ! write(*,*) 'co2: count=',count231 endif232 ! if (noms(iq).eq."co2_ice") then233 ! igcm_co2_ice=iq234 ! mmol(igcm_co2_ice)=44.235 ! count=count+1236 ! ! write(*,*) 'co2_ice: count=',count237 ! endif238 if (noms(iq).eq."h2o_vap") then239 igcm_h2o_vap=iq240 mmol(igcm_h2o_vap)=18.241 count=count+1242 ! write(*,*) 'h2o_vap: count=',count243 endif244 if (noms(iq).eq."h2o_ice") then245 igcm_h2o_ice=iq246 mmol(igcm_h2o_ice)=18.247 count=count+1248 ! write(*,*) 'h2o_ice: count=',count249 endif250 if (noms(iq).eq."co") then251 igcm_co=iq252 mmol(igcm_co)=28.253 count=count+1254 endif255 if (noms(iq).eq."o") then256 igcm_o=iq257 mmol(igcm_o)=16.258 count=count+1259 endif260 if (noms(iq).eq."o1d") then261 igcm_o1d=iq262 mmol(igcm_o1d)=16.263 count=count+1264 endif265 if (noms(iq).eq."o2") then266 igcm_o2=iq267 mmol(igcm_o2)=32.268 count=count+1269 endif270 if (noms(iq).eq."o3") then271 igcm_o3=iq272 mmol(igcm_o3)=48.273 count=count+1274 endif275 if (noms(iq).eq."h") then276 igcm_h=iq277 mmol(igcm_h)=1.278 count=count+1279 endif280 if (noms(iq).eq."h2") then281 igcm_h2=iq282 mmol(igcm_h2)=2.283 count=count+1284 endif285 if (noms(iq).eq."oh") then286 igcm_oh=iq287 mmol(igcm_oh)=17.288 count=count+1289 endif290 if (noms(iq).eq."ho2") then291 igcm_ho2=iq292 mmol(igcm_ho2)=33.293 count=count+1294 endif295 if (noms(iq).eq."h2o2") then296 igcm_h2o2=iq297 mmol(igcm_h2o2)=34.298 count=count+1299 endif300 225 if (noms(iq).eq."n2") then 301 226 igcm_n2=iq 302 227 mmol(igcm_n2)=28. 303 228 count=count+1 304 endif 305 if (noms(iq).eq."n2_ice") then 306 igcm_n2_ice=iq 307 mmol(igcm_n2_ice)=28. 308 count=count+1 309 ! write(*,*) 'n2_ice: count=',count 310 endif 311 if (noms(iq).eq."ch4") then 312 igcm_ch4=iq 313 mmol(igcm_ch4)=16. 314 count=count+1 315 endif 316 if (noms(iq).eq."ar") then 317 igcm_ar=iq 318 mmol(igcm_ar)=40. 319 count=count+1 320 endif 321 if (noms(iq).eq."n") then 322 igcm_n=iq 323 mmol(igcm_n)=14. 324 count=count+1 325 endif 326 if (noms(iq).eq."no") then 327 igcm_no=iq 328 mmol(igcm_no)=30. 329 count=count+1 330 endif 331 if (noms(iq).eq."no2") then 332 igcm_no2=iq 333 mmol(igcm_no2)=46. 334 count=count+1 335 endif 336 if (noms(iq).eq."n2d") then 337 igcm_n2d=iq 338 mmol(igcm_n2d)=28. 339 count=count+1 340 endif 341 if (noms(iq).eq."ch3") then 342 igcm_ch3=iq 343 mmol(igcm_ch3)=15. 344 count=count+1 345 endif 346 if (noms(iq).eq."ch") then 347 igcm_ch=iq 348 mmol(igcm_ch)=13. 349 count=count+1 350 endif 351 if (noms(iq).eq."3ch2") then 352 igcm_3ch2=iq 353 mmol(igcm_3ch2)=14. 354 count=count+1 355 endif 356 if (noms(iq).eq."1ch2") then 357 igcm_1ch2=iq 358 mmol(igcm_1ch2)=14. 359 count=count+1 360 endif 361 if (noms(iq).eq."cho") then 362 igcm_cho=iq 363 mmol(igcm_cho)=29. 364 count=count+1 365 endif 366 if (noms(iq).eq."ch2o") then 367 igcm_ch2o=iq 368 mmol(igcm_ch2o)=30. 369 count=count+1 370 endif 371 if (noms(iq).eq."ch3o") then 372 igcm_ch3o=iq 373 mmol(igcm_ch3o)=31. 374 count=count+1 375 endif 376 if (noms(iq).eq."c") then 377 igcm_c=iq 378 mmol(igcm_c)=12. 379 count=count+1 380 endif 381 if (noms(iq).eq."c2") then 382 igcm_c2=iq 383 mmol(igcm_c2)=24. 384 count=count+1 385 endif 386 if (noms(iq).eq."c2h") then 387 igcm_c2h=iq 388 mmol(igcm_c2h)=25. 389 count=count+1 390 endif 391 if (noms(iq).eq."c2h2") then 392 igcm_c2h2=iq 393 mmol(igcm_c2h2)=26. 394 count=count+1 395 endif 396 if (noms(iq).eq."c2h3") then 397 igcm_c2h3=iq 398 mmol(igcm_c2h3)=27. 399 count=count+1 400 endif 401 if (noms(iq).eq."c2h4") then 402 igcm_c2h4=iq 403 mmol(igcm_c2h4)=28. 404 count=count+1 405 endif 406 if (noms(iq).eq."c2h6") then 407 igcm_c2h6=iq 408 mmol(igcm_c2h6)=30. 409 count=count+1 410 endif 411 if (noms(iq).eq."ch2co") then 412 igcm_ch2co=iq 413 mmol(igcm_ch2co)=42. 414 count=count+1 415 endif 416 if (noms(iq).eq."ch3co") then 417 igcm_ch3co=iq 418 mmol(igcm_ch3co)=43. 419 count=count+1 420 endif 421 if (noms(iq).eq."hcaer") then 422 igcm_hcaer=iq 423 mmol(igcm_hcaer)=50. 424 count=count+1 425 endif 426 229 write(*,*) 'Tracer ',count,' = n2' 230 endif 231 if (noms(iq).eq."ch4_gas") then 232 igcm_ch4_gas=iq 233 mmol(igcm_ch4_gas)=16. 234 count=count+1 235 write(*,*) 'Tracer ',count,' = ch4 gas' 236 endif 237 if (noms(iq).eq."ch4_ice") then 238 igcm_ch4_ice=iq 239 mmol(igcm_ch4_ice)=16. 240 radius(iq)=3.e-6 241 rho_q(iq)=rho_ch4_ice 242 count=count+1 243 write(*,*) 'Tracer ',count,' = ch4 ice' 244 endif 245 if (noms(iq).eq."co_gas") then 246 igcm_co_gas=iq 247 mmol(igcm_co_gas)=28. 248 count=count+1 249 write(*,*) 'Tracer ',count,' = co gas' 250 endif 251 if (noms(iq).eq."co_ice") then 252 igcm_co_ice=iq 253 mmol(igcm_co_ice)=28. 254 radius(iq)=3.e-6 255 rho_q(iq)=rho_co_ice 256 count=count+1 257 write(*,*) 'Tracer ',count,' = co ice' 258 endif 259 if (noms(iq).eq."prec_haze") then 260 igcm_prec_haze=iq 261 count=count+1 262 write(*,*) 'Tracer ',count,' = prec haze' 263 endif 264 if (noms(iq).eq."haze") then 265 igcm_haze=iq 266 count=count+1 267 radius(iq)=rad_haze 268 rho_q(iq)=rho_haze 269 write(*,*) 'Tracer ',count,' = haze' 270 endif 271 if (noms(iq).eq."haze10") then 272 igcm_haze10=iq 273 count=count+1 274 radius(iq)=10.e-9 275 rho_q(iq)=rho_haze 276 write(*,*) 'Tracer ',count,' = haze10' 277 endif 278 if (noms(iq).eq."haze30") then 279 igcm_haze30=iq 280 count=count+1 281 radius(iq)=30.e-9 282 rho_q(iq)=rho_haze 283 write(*,*) 'Tracer ',count,' = haze30' 284 endif 285 if (noms(iq).eq."haze50") then 286 igcm_haze50=iq 287 count=count+1 288 radius(iq)=50.e-9 289 rho_q(iq)=rho_haze 290 write(*,*) 'Tracer ',count,' = haze50' 291 endif 292 if (noms(iq).eq."haze100") then 293 igcm_haze100=iq 294 count=count+1 295 radius(iq)=100.e-9 296 rho_q(iq)=rho_haze 297 write(*,*) 'Tracer ',count,' = haze100' 298 endif 299 ! Eddy diffusion tracers 300 if (noms(iq).eq."eddy1e6") then 301 igcm_eddy1e6=iq 302 count=count+1 303 write(*,*) 'Tracer ',count,' = eddy1e6' 304 endif 305 if (noms(iq).eq."eddy1e7") then 306 igcm_eddy1e7=iq 307 count=count+1 308 write(*,*) 'Tracer ',count,' = eddy1e7' 309 endif 310 if (noms(iq).eq."eddy5e7") then 311 igcm_eddy5e7=iq 312 count=count+1 313 write(*,*) 'Tracer ',count,' = eddy5e7' 314 endif 315 if (noms(iq).eq."eddy1e8") then 316 igcm_eddy1e8=iq 317 count=count+1 318 write(*,*) 'Tracer ',count,' = eddy1e8' 319 endif 320 if (noms(iq).eq."eddy5e8") then 321 igcm_eddy5e8=iq 322 count=count+1 323 write(*,*) 'Tracer ',count,' = eddy5e8' 324 endif 427 325 enddo ! of do iq=1,nq 428 326 429 327 ! 3. find condensable traceurs different from h2o and n2 430 328 do iq=1,nq 431 if ((index(noms(iq),"vap") .ne. 0) .and. (index(noms(iq)," h2o") .eq. 0) .and. (index(noms(iq),"n2") .eq. 0)) then432 count=count+1 433 endif 434 if ((index(noms(iq),"ice") .ne. 0) .and. (index(noms(iq)," h2o") .eq. 0) .and. (index(noms(iq),"n2") .eq. 0)) then329 if ((index(noms(iq),"vap") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then 330 count=count+1 331 endif 332 if ((index(noms(iq),"ice") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then 435 333 count=count+1 436 334 endif 437 335 438 336 enddo ! of do iq=1,nq 439 337 440 338 ! check that we identified all tracers: 441 339 if (count.ne.nq) then … … 456 354 if (is_master) then 457 355 rewind(407) 458 do 356 do 459 357 read(407,'(A)') tracline 460 if (index(tracline,"#") .ne. 1) then 358 if (index(tracline,"#") .ne. 1) then 461 359 exit 462 360 endif … … 502 400 n_rgcs = sum(is_rgcs) 503 401 write(*,*)'Number of Radiative Generic Condensable Species is n_rgcs = ',n_rgcs 504 if (n_rgcs> ngt/2) then 402 if (n_rgcs> ngt/2) then 505 403 write(*,*) 'You have more Radiative Generic Condensable Species than Generic Condensable Species' 506 404 write(*,*)'This is not possible: check your Modern traceur.def' … … 523 421 lw_ch4=586700. 524 422 lw_n2=2.5e5 525 423 526 424 if (haze) then 527 425 ! the sedimentation radius remains radius(igcm_haze) … … 530 428 else 531 429 nmono=1 532 endif 430 endif 533 431 534 432 ia=0 … … 539 437 540 438 block=0 ! Only one type of haze is active : the first one set in traceur.def 541 do iq=1,nq mx439 do iq=1,nq 542 440 tracername=noms(iq) 543 write(*,*) "--> tracername ",iq,'/',nq mx,' = ',tracername441 write(*,*) "--> tracername ",iq,'/',nq,' = ',tracername 544 442 if (tracername(1:4).eq."haze".and.block.eq.0) then 545 443 i_haze=iq 546 444 block=1 547 445 write(*,*) "i_haze=",i_haze 548 write(*,*) "Careful: if you set many haze traceurs in 549 & traceur.def,only ",tracername," will be radiatively active550 &(first one in traceur.def)"446 write(*,*) "Careful: if you set many haze traceurs in& 447 traceur.def,only ",tracername," will be radiatively active& 448 (first one in traceur.def)" 551 449 endif 552 450 enddo … … 563 461 write(*,*) 'alpha_devil = ', alpha_devil 564 462 write(*,*) 'radius = ', radius 565 write(*,*) 'Qext = ', qext 463 write(*,*) 'Qext = ', qext 566 464 write(*,*) 567 465 … … 579 477 integer, intent(in) :: iq ! tracer index 580 478 character(len=500),intent(in) :: tracline ! traceur.def lines with parameters 581 479 582 480 read(tracline,*) noms(iq) 583 481 ! JVO 20 : We should add a sanity check aborting when duplicates in names ! … … 593 491 mmol(iq) 594 492 end if 595 ! option aki 493 ! option aki 596 494 if (index(tracline,'aki=' ) /= 0) then 597 495 read(tracline(index(tracline,'aki=')+len('aki='):),*) & … … 603 501 aki(iq) 604 502 end if 605 ! option cpi 503 ! option cpi 606 504 if (index(tracline,'cpi=' ) /= 0) then 607 505 read(tracline(index(tracline,'cpi=')+len('cpi='):),*) & … … 613 511 cpi(iq) 614 512 end if 615 ! option is_chim 513 ! option is_chim 616 514 if (index(tracline,'is_chim=' ) /= 0) then 617 515 read(tracline(index(tracline,'is_chim=')+len('is_chim='):),*) & … … 623 521 is_chim(iq) 624 522 end if 625 ! option is_rad 523 ! option is_rad 626 524 if (index(tracline,'is_rad=') /= 0) then 627 525 read(tracline(index(tracline,'is_rad=') & … … 668 566 end if 669 567 !option is_condensable (LT) 670 if (index(tracline,'is_condensable=') /=0) then 568 if (index(tracline,'is_condensable=') /=0) then 671 569 read(tracline(index(tracline,'is_condensable=') & 672 570 +len('is_condensable='):),*) is_condensable(iq) 673 571 write(*,*) ' Parameter value (traceur.def) :'// & 674 ' is_condensable=', is_condensable(iq) 572 ' is_condensable=', is_condensable(iq) 675 573 else 676 574 write(*,*) ' Parameter value (default) :'// & 677 ' is_condensable=', is_condensable(iq) 575 ' is_condensable=', is_condensable(iq) 678 576 endif 679 577 !option radius 680 if (index(tracline,'radius=') .ne. 0) then 578 if (index(tracline,'radius=') .ne. 0) then 681 579 read(tracline(index(tracline,'radius=') & 682 580 +len('radius='):),*) radius(iq) … … 685 583 else 686 584 write(*,*) ' Parameter value (default) :'// & 687 ' radius=', radius(iq)," m" 585 ' radius=', radius(iq)," m" 688 586 endif 689 587 !option rho 690 if (index(tracline,'rho=') .ne. 0) then 588 if (index(tracline,'rho=') .ne. 0) then 691 589 read(tracline(index(tracline,'rho=') & 692 590 +len('rho='):),*) rho_q(iq) … … 695 593 else 696 594 write(*,*) ' Parameter value (default) :'// & 697 ' rho=', rho_q(iq) 595 ' rho=', rho_q(iq) 698 596 endif 699 !option is_rgcs 700 if (index(tracline,'is_rgcs') .ne. 0) then 597 !option is_rgcs 598 if (index(tracline,'is_rgcs') .ne. 0) then 701 599 read(tracline(index(tracline,'is_rgcs=') & 702 600 +len('is_rgcs='):),*) is_rgcs(iq) -
trunk/LMDZ.PLUTO/libf/phypluto/phys_state_var_mod.F90
r3184 r3195 13 13 use comsaison_h, only: mu0, fract 14 14 use radcommon_h, only: glat 15 ! use slab_ice_h, only : noceanmx16 ! USE ocean_slab_mod, ONLY: nslay17 15 USE radinc_h, only : L_NSPECTI, L_NSPECTV,naerkind 18 16 use surfdat_h, only: phisfi, albedodat, & … … 73 71 real,dimension(:,:,:),allocatable,save :: int_dtauv ! VI optical thickness of layers within narrowbands for diags (). 74 72 real,dimension(:,:,:),allocatable,save :: int_dtaui ! IR optical thickness of layers within narrowbands for diags (). 75 !$OMP THREADPRIVATE(int_dtaui,int_dtauv) 73 !$OMP THREADPRIVATE(int_dtaui,int_dtauv) 76 74 77 75 real,allocatable,dimension(:),save :: tau_col ! Total Aerosol Optical Depth. … … 100 98 real,dimension(:,:),allocatable,save :: lscaledEz ! Heat from largescale (W.m-2) 101 99 !$OMP THREADPRIVATE(madjdE,madjdEz,lscaledE,lscaledEz) 102 100 103 101 real,dimension(:),allocatable,save :: genericconddE ! Heat from generic condensation (W.m-2) 104 102 !$OMP THREADPRIVATE(genericconddE) … … 122 120 !zpic(:) ! Maximum de l'OESM 123 121 !zval(:) ! Minimum de l'OESM 124 !rugoro(:) ! longueur de rugosite de l'OESM 122 !rugoro(:) ! longueur de rugosite de l'OESM 125 123 126 124 ALLOCATE(phisfi(klon)) -
trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
r3184 r3195 1 1 module physiq_mod 2 2 3 3 implicit none 4 4 5 5 contains 6 6 7 7 subroutine physiq(ngrid,nlayer,nq, & 8 8 firstcall,lastcall, & … … 14 14 15 15 !! 16 use write_field_phy, only: Writefield_phy 16 use write_field_phy, only: Writefield_phy 17 17 !! 18 use ioipsl_getin_p_mod, only: getin_p 18 use ioipsl_getin_p_mod, only: getin_p 19 19 use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir 20 20 use generic_cloud_common_h, only : epsi_generic, Psat_generic … … 23 23 use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV 24 24 use suaer_corrk_mod, only: suaer_corrk 25 use radii_mod, only: n2_reffrad26 use aerosol_mod, only: i niaerosol, iaero_n225 use radii_mod, only: su_aer_radii,haze_reffrad_fix 26 use aerosol_mod, only: iaero_haze, i_haze, haze_prof 27 27 use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, & 28 28 dryness … … 33 33 USE comgeomfi_h, only: totarea, totarea_planet 34 34 USE tracer_h, only: noms, mmol, radius, rho_q, qext, & 35 igcm_n2,& 35 36 alpha_lift, alpha_devil, qextrhor, & 36 igcm_dustbin, & 37 igcm_n2_ice, nesp, is_chim, is_condensable,constants_epsi_generic 37 nesp, is_chim, is_condensable,constants_epsi_generic 38 38 use time_phylmdz_mod, only: ecritphy, iphysiq, nday 39 39 use phyetat0_mod, only: phyetat0 … … 48 48 use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, & 49 49 callrad, callsoil, nosurf, & 50 n2cond, corrk, diagdtau,&50 aerohaze, corrk, diagdtau,& 51 51 diurnal, enertest, fat1au, & 52 52 icetstep, intheat, iradia, kastprof, & 53 53 lwrite, mass_redistrib, meanOLR, & 54 n earn2cond, noseason_day, &54 n2cond,nearn2cond, noseason_day, & 55 55 season, sedimentation,generic_condensation, & 56 aerohaze, haze_proffix, & 56 57 specOLR, & 57 58 startphy_file, testradtimes, & … … 67 68 use turb_mod, only : q2,sensibFlux,turb_resolved 68 69 use mass_redistribution_mod, only: mass_redistribution 69 use condensation_generic_mod, only: condensation_generic 70 use condensation_generic_mod, only: condensation_generic 70 71 use datafile_mod, only: datadir 71 72 #ifndef MESOSCALE … … 85 86 #endif 86 87 87 #ifdef CPP_XIOS 88 #ifdef CPP_XIOS 88 89 use xios_output_mod, only: initialize_xios_output, & 89 90 update_xios_timestep, & … … 96 97 97 98 !================================================================== 98 ! 99 ! 99 100 ! Purpose 100 101 ! ------- … … 186 187 ! Frederic Hourdin 15/10/93 187 188 ! Francois Forget 1994 188 ! Christophe Hourdin 02/1997 189 ! Christophe Hourdin 02/1997 189 190 ! Subroutine completely rewritten by F. Forget (01/2000) 190 191 ! Water ice clouds: Franck Montmessin (update 06/2003) 191 192 ! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009) 192 193 ! New correlated-k radiative scheme: R. Wordsworth (2009) 193 ! Many specifically Martian subroutines removed: R. Wordsworth (2009) 194 ! Many specifically Martian subroutines removed: R. Wordsworth (2009) 194 195 ! Improved water cycle: R. Wordsworth / B. Charnay (2010) 195 196 ! To F90: R. Wordsworth (2010) 196 197 ! New turbulent diffusion scheme: J. Leconte (2012) 197 198 ! Loops converted to F90 matrix format: J. Leconte (2012) 198 ! No more ngrid mx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)199 ! No more ngrid/nq, F90 commons and adaptation to parallel: A. Spiga (2012) 199 200 ! Purge of the code : M. Turbet (2015) 200 201 ! Photochemical core developped by F. Lefevre: B. Charnay (2017) … … 217 218 integer,intent(in) :: nlayer ! Number of atmospheric layers. 218 219 integer,intent(in) :: nq ! Number of tracers. 219 220 220 221 logical,intent(in) :: firstcall ! Signals first call to physics. 221 222 logical,intent(in) :: lastcall ! Signals last call to physics. 222 223 223 224 real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0. 224 225 real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon). … … 253 254 ! ----------------- 254 255 255 ! Aerosol (dust or ice) extinction optical depth at reference wavelength 256 ! Aerosol (dust or ice) extinction optical depth at reference wavelength 256 257 ! for the "naerkind" optically active aerosols: 257 258 … … 264 265 integer l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice 265 266 logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer 266 267 267 268 real zls ! Solar longitude (radians). 268 269 real zlss ! Sub solar point longitude (radians). … … 272 273 273 274 ! VARIABLES for the thermal plume model (AF24: deleted) 274 275 275 276 ! TENDENCIES due to various processes : 276 277 … … 281 282 real zdtsdif(ngrid) ! Turbdiff/vdifc routines. 282 283 ! real zdtsurf_hyd(ngrid) ! Hydrol routine. 283 284 ! For Atmospheric Temperatures : (K/s) 284 285 ! For Atmospheric Temperatures : (K/s) 285 286 real dtlscale(ngrid,nlayer) ! Largescale routine. 286 real dt_generic_condensation(ngrid,nlayer) ! condensation_generic routine. 287 real dt_generic_condensation(ngrid,nlayer) ! condensation_generic routine. 287 288 real zdtc(ngrid,nlayer) ! Condense_n2 routine. 288 289 real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines. … … 290 291 real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine. 291 292 real zdtchim(ngrid,nlayer) ! Calchim routine. 292 293 293 294 ! For Surface Tracers : (kg/m2/s) 294 295 real dqsurf(ngrid,nq) ! Cumulated tendencies. 295 296 real zdqsurfc(ngrid) ! Condense_n2 routine. 297 REAL zdqsc(ngrid,nq) ! Condense_n2 routine. 296 298 real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines. 297 299 real zdqssed(ngrid,nq) ! Callsedim routine. 298 300 real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine. 299 301 300 302 ! For Tracers : (kg/kg_of_air/s) 301 303 real zdqc(ngrid,nlayer,nq) ! Condense_n2 routine. … … 311 313 !$OMP THREADPRIVATE(zdqchim,zdqschim) 312 314 315 316 !! PLUTO variables 317 REAL zdqch4cloud(ngrid,nlayer,nq) 318 REAL zdqsch4cloud(ngrid,nq) 319 REAL zdtch4cloud(ngrid,nlayer) 320 REAL zdqcocloud(ngrid,nlayer,nq) 321 REAL zdqscocloud(ngrid,nq) 322 REAL zdtcocloud(ngrid,nlayer) 323 REAL rice_ch4(ngrid,nlayer) ! Methane ice geometric mean radius (m) 324 REAL rice_co(ngrid,nlayer) ! CO ice geometric mean radius (m) 325 326 REAL zdqsch4fast(ngrid) ! used only for fast model nogcm 327 REAL zdqch4fast(ngrid) ! used only for fast model nogcm 328 REAL zdqscofast(ngrid) ! used only for fast model nogcm 329 REAL zdqcofast(ngrid) ! used only for fast model nogcm 330 REAL zdqflow(ngrid,nq) 331 332 REAL zdteuv(ngrid,nlayer) ! (K/s) 333 REAL zdtconduc(ngrid,nlayer) ! (K/s) 334 REAL zdumolvis(ngrid,nlayer) 335 REAL zdvmolvis(ngrid,nlayer) 336 real zdqmoldiff(ngrid,nlayer,nq) 337 338 ! Haze relatated tendancies 339 REAL zdqhaze(ngrid,nlayer,nq) 340 REAL zdqprodhaze(ngrid,nq) 341 REAL zdqsprodhaze(ngrid) 342 REAL zdqhaze_col(ngrid) 343 REAL zdqphot_prec(ngrid,nlayer) 344 REAL zdqphot_ch4(ngrid,nlayer) 345 REAL zdqconv_prec(ngrid,nlayer) 346 REAL zdq_source(ngrid,nlayer,nq) 347 ! Fast Haze relatated tendancies 348 REAL fluxbot(ngrid) 349 REAL gradflux(ngrid) 350 REAL fluxlym_sol_bot(ngrid) ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface 351 REAL fluxlym_ipm_bot(ngrid) ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface 352 REAL flym_sol(ngrid) ! Incident Solar flux Lyman alpha ph.m-2.s-1 353 REAL flym_ipm(ngrid) ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 354 355 356 313 357 REAL array_zero1(ngrid) 314 358 REAL array_zero2(ngrid,nlayer) 315 359 316 360 ! For Winds : (m/s/s) 317 361 real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer) ! Convadj routine. … … 320 364 real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines. 321 365 real zdhadj(ngrid,nlayer) ! Convadj routine. 322 366 REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer) ! condense_n2 routine. 367 323 368 ! For Pressure and Mass : 324 369 real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s). … … 326 371 real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s). 327 372 328 329 373 374 330 375 ! Local variables for LOCAL CALCULATIONS: 331 376 ! --------------------------------------- … … 352 397 real vmr(ngrid,nlayer) ! volume mixing ratio 353 398 real time_phys 354 399 355 400 real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic. 356 401 357 402 real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2). 403 404 ! Pluto adding variables 405 real vmr_ch4(ngrid) ! vmr ch4 406 real vmr_co(ngrid) ! vmr co 407 real rho(ngrid,nlayer) ! density 408 real zrho_ch4(ngrid,nlayer) ! density methane kg.m-3 409 real zrho_co(ngrid,nlayer) ! density CO kg.m-3 410 real zrho_haze(ngrid,nlayer) ! density haze kg.m-3 411 real zdqrho_photprec(ngrid,nlayer) !photolysis rate kg.m-3.s-1 412 real zq1temp_ch4(ngrid) ! 413 real qsat_ch4(ngrid) ! 414 real qsat_ch4_l1(ngrid) ! 415 ! CHARACTER(LEN=20) :: txt ! to temporarly store text for eddy tracers 416 real profmmr(ngrid,nlayer) ! fixed profile of haze if haze_proffix 417 real sensiblehf1(ngrid) ! sensible heat flux 418 real sensiblehf2(ngrid) ! sensible heat flux 358 419 359 420 ! included by RW for H2O Manabe scheme … … 367 428 real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW 368 429 !$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW) 369 430 370 431 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily 371 432 372 real dtmoist_max,dtmoist_min 433 real dtmoist_max,dtmoist_min 373 434 real dItot, dItot_tmp, dVtot, dVtot_tmp 374 435 … … 408 469 real :: flux_sens_lat(ngrid) 409 470 real :: qsurfint(ngrid,nq) 410 #ifdef MESOSCALE 471 #ifdef MESOSCALE 411 472 REAL :: lsf_dt(nlayer) 412 473 REAL :: lsf_dq(nlayer) … … 416 477 logical, save :: check_physics_inputs=.false. 417 478 logical, save :: check_physics_outputs=.false. 418 !$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs) 479 !$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs) 419 480 420 481 ! Misc 421 482 character*2 :: str2 422 character(len=10) :: tmp1 483 character(len=10) :: tmp1 423 484 character(len=10) :: tmp2 424 485 !================================================================================================== … … 458 519 ! Initialize aerosol indexes. 459 520 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 460 call iniaerosol521 ! call iniaerosol 461 522 ! allocate related local arrays 462 523 ! (need be allocated instead of automatic because of "naerkind") … … 472 533 ! Read 'startfi.nc' file. 473 534 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 474 #ifndef MESOSCALE 535 #ifndef MESOSCALE 475 536 call phyetat0(startphy_file, & 476 537 ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, & 477 538 day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf) 478 479 #else 539 540 #else 480 541 481 542 day_ini = pday … … 505 566 write (*,*) 'In physiq day_ini =', day_ini 506 567 507 ! Initialize albedo calculation. 568 ! Initialize albedo calculation. 508 569 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 509 570 albedo(:,:)=0.0 … … 511 572 albedo_snow_SPECTV(:)=0.0 512 573 albedo_n2_ice_SPECTV(:)=0.0 513 call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV) 514 515 ! Initialize orbital calculation. 574 call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV) 575 576 ! Initialize orbital calculation. 516 577 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 517 578 call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) … … 555 616 endif 556 617 557 if(meanOLR)then 558 call system('rm -f rad_bal.out') ! to record global radiative balance. 559 call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures. 618 if(meanOLR)then 619 call system('rm -f rad_bal.out') ! to record global radiative balance. 620 call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures. 560 621 call system('rm -f h2o_bal.out') ! to record global hydrological balance. 561 622 endif … … 563 624 564 625 !! Initialize variables for water cycle ! AF24: removed 565 626 566 627 ! Set metallicity for GCS 567 628 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 568 629 metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic 569 630 call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic 570 631 571 632 ! Set some parameters for the thermal plume model !AF24: removed 572 633 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 573 634 574 635 #ifndef MESOSCALE 575 636 if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d. … … 580 641 581 642 !! call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1) 582 #endif 583 if (corrk) then 584 ! We initialise the spectral grid here instead of 585 ! at firstcall of callcorrk so we can output XspecIR, XspecVI 643 #endif 644 if (corrk) then 645 ! We initialise the spectral grid here instead of 646 ! at firstcall of callcorrk so we can output XspecIR, XspecVI 586 647 ! when using Dynamico 587 648 print*, "physiq_mod: Correlated-k data base folder:",trim(datadir) … … 595 656 call setspv ! Basic visible properties. 596 657 call sugas_corrk ! Set up gaseous absorption properties. 597 call suaer_corrk ! Set up aerosol optical properties. 658 if (aerohaze) then 659 call suaer_corrk ! Set up aerosol optical properties. 660 endif 598 661 endif 599 662 … … 606 669 year_day,presnivs,pseudoalt,WNOI,WNOV) 607 670 #endif 608 671 609 672 !! call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1) 610 673 … … 613 676 614 677 !! call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1) 615 !! call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap)) 678 !! call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap)) 616 679 617 680 if (check_physics_inputs) then … … 627 690 ! ------------------------------------------------------ 628 691 629 #ifdef CPP_XIOS 692 #ifdef CPP_XIOS 630 693 ! update XIOS time/calendar 631 694 call update_xios_timestep 632 #endif 695 #endif 633 696 634 697 ! Initialize various variables 635 698 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 636 699 637 700 if ( .not.nearn2cond ) then 638 701 pdt(1:ngrid,1:nlayer) = 0.0 639 endif 702 endif 640 703 zdtsurf(1:ngrid) = 0.0 641 704 pdq(1:ngrid,1:nlayer,1:nq) = 0.0 … … 644 707 pdv(1:ngrid,1:nlayer) = 0.0 645 708 pdpsrf(1:ngrid) = 0.0 646 zflubid(1:ngrid) = 0.0 709 zflubid(1:ngrid) = 0.0 647 710 flux_sens_lat(1:ngrid) = 0.0 648 711 taux(1:ngrid) = 0.0 … … 660 723 661 724 call orbite(zls,dist_star,declin,right_ascen) 662 725 663 726 if (diurnal) then 664 727 zlss=-2.*pi*(zday-.5) 665 728 else if(diurnal .eqv. .false.) then 666 729 zlss=9999. 667 endif 730 endif 668 731 669 732 … … 673 736 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 674 737 zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer) 675 do l=1,nlayer 738 do l=1,nlayer 676 739 zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid) 677 740 enddo … … 685 748 zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) 686 749 enddo 687 enddo 750 enddo 688 751 689 752 !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22 690 753 691 754 zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1) 692 755 … … 694 757 ! Note : Potential temperature calculation may not be the same in physiq and dynamic... 695 758 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 696 do l=1,nlayer 759 do l=1,nlayer 697 760 do ig=1,ngrid 698 761 zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp … … 745 808 endif 746 809 747 endif 810 endif 748 811 749 812 if (callrad) then 750 813 if( mod(icount-1,iradia).eq.0.or.lastcall) then 751 814 752 815 ! Eclipse incoming sunlight !AF24: removed 753 816 754 817 !! call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap)) 755 !! call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap)) 818 !! call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap)) 756 819 757 820 758 821 if (corrk) then 759 822 760 823 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 761 824 ! II.a Call correlated-k radiative transfer scheme … … 766 829 endif 767 830 ! if(water) then !AF24: removed 768 831 769 832 if(generic_condensation) then 770 833 do iq=1,nq 771 834 772 835 call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic) 773 836 774 837 if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer 775 838 776 839 epsi_generic=constants_epsi_generic(iq) 777 840 778 muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap)) 779 muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap)) 841 muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap)) 842 muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap)) 780 843 781 844 endif 782 end do ! do iq=1,nq loop on tracers 845 end do ! do iq=1,nq loop on tracers 783 846 ! take into account generic condensable specie (GCS) effect on mean molecular weight 784 847 785 848 else 786 849 muvar(1:ngrid,1:nlayer+1)=mugaz 787 endif 788 850 endif 851 789 852 ! if(ok_slab_ocean) then !AF24: removed 790 853 791 854 ! standard callcorrk 792 855 ! clearsky=.false. … … 799 862 int_dtaui,int_dtauv, & 800 863 tau_col,cloudfrac,totcloudfrac, & 801 .false.,firstcall,lastcall) 864 .false.,firstcall,lastcall) 802 865 803 866 !GG (feb2021): Option to "artificially" decrease the raditive time scale in 804 867 !the deep atmosphere press > 0.1 bar. Suggested by MT. 805 !! COEFF_RAD to be "tuned" to facilitate convergence of tendency 806 868 !! COEFF_RAD to be "tuned" to facilitate convergence of tendency 869 807 870 !coeff_rad=0. ! 0 values, it doesn't accelerate the convergence 808 871 !coeff_rad=0.5 809 !coeff_rad=1. 872 !coeff_rad=1. 810 873 !do l=1, nlayer 811 874 ! do ig=1,ngrid … … 836 899 ! Net atmospheric radiative heating rate (K.s-1) 837 900 dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer) 838 839 ! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that ! 901 902 ! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that ! 840 903 if (firstcall .and. albedo_spectral_mode) then 841 904 call spectral_albedo_calc(albedo_snow_SPECTV) … … 847 910 else 848 911 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 849 ! II.c Atmosphere has no radiative effect 912 ! II.c Atmosphere has no radiative effect 850 913 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 851 914 fluxtop_dn(1:ngrid) = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2 … … 856 919 print*,'------------WARNING---WARNING------------' ! by MT2015. 857 920 print*,'You are in corrk=false mode, ' 858 print*,'and the surface albedo is taken equal to the first visible spectral value' 859 921 print*,'and the surface albedo is taken equal to the first visible spectral value' 922 860 923 fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1)) 861 924 fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid) … … 867 930 868 931 endif ! of if(mod(icount-1,iradia).eq.0) 869 932 870 933 871 934 ! Transformation of the radiative tendencies … … 875 938 fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid) 876 939 pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer) 877 940 878 941 ! Test of energy conservation 879 942 !---------------------------- … … 908 971 909 972 if (calldifv) then 910 973 911 974 zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid) 912 975 913 976 ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff. 914 977 if (UseTurbDiff) then 915 978 916 979 call turbdiff(ngrid,nlayer,nq, & 917 980 ptimestep,capcal, & … … 924 987 925 988 else 926 989 927 990 zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer) 928 991 929 992 call vdifc(ngrid,nlayer,nq,zpopsk, & 930 993 ptimestep,capcal,lwrite, & … … 955 1018 !! call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_vap)) 956 1019 957 if (tracer) then 1020 if (tracer) then 958 1021 pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq) 959 1022 dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq) … … 966 1029 !------------------------- 967 1030 if(enertest)then 968 1031 969 1032 dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:) 970 1033 do ig = 1, ngrid … … 972 1035 dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground 973 1036 enddo 974 1037 975 1038 call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot) 976 1039 dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:) 977 1040 call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots) 978 1041 call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux) 979 1042 980 1043 if (is_master) then 981 1044 982 1045 if (UseTurbDiff) then 983 1046 print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2' … … 990 1053 end if 991 1054 endif ! end of 'is_master' 992 1055 993 1056 ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere. 994 1057 endif ! end of 'enertest' … … 1008 1071 ! IV. Convection : 1009 1072 !------------------- 1010 1073 1011 1074 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~ 1012 1075 ! IV.a Thermal plume model : AF24: removed 1013 1076 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~ 1014 1077 1015 1078 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1016 1079 ! IV.b Dry convective adjustment : … … 1037 1100 zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only 1038 1101 1039 if(tracer) then 1040 pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq) 1102 if(tracer) then 1103 pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq) 1041 1104 end if 1042 1105 … … 1048 1111 1049 1112 ! ! Test water conservation !AF24: removed 1050 1113 1051 1114 endif ! end of 'calladj' 1052 1115 !---------------------------------------------- 1053 1116 ! Non orographic Gravity Waves: AF24: removed 1054 1117 !--------------------------------------------- 1055 1118 1056 1119 !----------------------------------------------- 1057 ! V. Carbon dioxidecondensation-sublimation :1120 ! V. Nitrogen condensation-sublimation : 1058 1121 !----------------------------------------------- 1059 1122 1060 1123 if (n2cond) then 1061 1124 1062 1125 if (.not.tracer) then 1063 1126 print*,'We need a N2 ice tracer to condense N2' 1064 1127 call abort 1065 1128 endif 1129 1066 1130 call condense_n2(ngrid,nlayer,nq,ptimestep, & 1067 capcal,pplay,pplev,tsurf,pt, 1068 p dt,zdtsurf,pq,pdq,&1069 qsurf ,zdqsurfc,albedo,emis,&1070 albedo_bareground,albedo_n2_ice_SPECTV,&1071 zd tc,zdtsurfc,pdpsrf,zdqc)1131 capcal,pplay,pplev,tsurf,pt, & 1132 pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, & 1133 qsurf(1,igcm_n2),albedo,emis, & 1134 zdtc,zdtsurfc,pdpsrf,zduc,zdvc, & 1135 zdqc,zdqsc(1,igcm_n2)) 1072 1136 1073 1137 pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer) … … 1075 1139 1076 1140 pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq) 1077 dqsurf(1:ngrid,igcm_n2_ice) = dqsurf(1:ngrid,igcm_n2_ice) + zdqsurfc(1:ngrid)1141 ! dqsurf(1:ngrid,igcm_n2_ice) = dqsurf(1:ngrid,igcm_n2_ice) + zdqsurfc(1:ngrid) 1078 1142 1079 1143 !! call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap)) … … 1094 1158 1095 1159 !--------------------------------------------- 1096 ! VI. Specific parameterizations for tracers 1160 ! VI. Specific parameterizations for tracers 1097 1161 !--------------------------------------------- 1098 1162 1099 1163 if (tracer) then 1100 1164 1101 1165 ! --------------------- 1102 1166 ! VI.1. Water and ice !AF24: removed … … 1110 1174 1111 1175 !Generic Condensation 1112 if (generic_condensation) then 1176 if (generic_condensation) then 1113 1177 call condensation_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay, & 1114 1178 pt,pq,pdt,pdq,dt_generic_condensation, & … … 1128 1192 end if 1129 1193 1130 ! if (.not. water) then 1194 ! if (.not. water) then 1131 1195 ! Compute GCS (Generic Condensable Specie) cloud fraction. For now we can not have both water cloud fraction and GCS cloud fraction 1132 1196 ! Water is the priority … … 1150 1214 enddo 1151 1215 endif 1152 end do ! do iq=1,nq loop on tracers 1216 end do ! do iq=1,nq loop on tracers 1153 1217 ! endif ! .not. water 1154 1218 … … 1164 1228 1165 1229 ! if(watertest)then !AF24: removed 1166 1230 1167 1231 call callsedim(ngrid,nlayer,ptimestep, & 1168 1232 pplev,zzlev,pt,pdt,pq,pdq, & … … 1170 1234 1171 1235 ! if(watertest)then !AF24: removed 1172 1236 1173 1237 ! Whether it falls as rain or snow depends only on the surface temperature 1174 1238 pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq) … … 1178 1242 1179 1243 ! ! Test water conservation !AF24: removed 1180 1244 1181 1245 end if ! end of 'sedimentation' 1182 1246 … … 1200 1264 zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer)) 1201 1265 enddo 1202 1266 1203 1267 ! call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr) 1204 1268 ! call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col) … … 1209 1273 pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, & 1210 1274 zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr) 1211 1275 1212 1276 pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq) 1213 1277 dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq) … … 1217 1281 pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid) 1218 1282 zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid) 1219 1283 1220 1284 endif 1221 1285 … … 1238 1302 qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq) 1239 1303 1240 ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water 1304 ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water 1241 1305 ! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain. 1242 1306 qsurf_hist(:,:) = qsurf(:,:) … … 1250 1314 1251 1315 !------------------------------------------------ 1252 ! VII. Surface and sub-surface soil temperature 1316 ! VII. Surface and sub-surface soil temperature 1253 1317 !------------------------------------------------ 1254 1318 … … 1256 1320 ! ! Increment surface temperature 1257 1321 ! if(ok_slab_ocean)then !AF24: removed 1258 1259 tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) 1322 1323 tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) 1260 1324 1261 1325 ! Compute soil temperatures and subsurface heat flux. 1262 1326 if (callsoil) then 1263 1327 call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, & 1264 ptimestep,tsurf,tsoil,capcal,fluxgrd) 1328 ptimestep,tsurf,tsoil,capcal,fluxgrd) 1265 1329 endif 1266 1330 1267 1331 1268 1332 ! if (ok_slab_ocean) then !AF24: removed 1269 1333 1270 1334 ! Test energy conservation 1271 1335 if(enertest)then 1272 call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots) 1336 call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots) 1273 1337 if (is_master) print*,'Surface energy change =',dEtots,' W m-2' 1274 1338 endif … … 1286 1350 zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep 1287 1351 zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep 1288 1352 1289 1353 !! Recast thermal plume vertical velocity array for outputs 1290 1354 !! AF24: removed 1291 1355 1292 1356 ! Diagnostic. 1293 1357 zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep … … 1352 1416 endif 1353 1417 end do 1354 1418 1355 1419 if(ngrid.eq.1)then 1356 1420 DYN=0.0 1357 1421 endif 1358 1422 1359 1423 if (is_master) then 1360 1424 print*,' ISR ASR OLR GND DYN [W m^-2]' … … 1412 1476 ! Generalised for arbitrary aerosols now. By LK 1413 1477 reffcol(1:ngrid,1:naerkind)=0.0 1414 if(n2cond.and.(iaero_n2.ne.0))then 1415 call n2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_n2)) 1416 do ig=1,ngrid 1417 reffcol(ig,iaero_n2) = SUM(zq(ig,1:nlayer,igcm_n2_ice)*reffrad(ig,1:nlayer,iaero_n2)*mass(ig,1:nlayer)) 1418 enddo 1419 endif 1420 ! if(water.and.(iaero_h2o.ne.0))then !AF24: removed 1421 1478 ! call n2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_haze)) 1479 if (haze_proffix.and.i_haze.gt.0.) then 1480 call haze_prof(ngrid,nlayer,zzlay,pplay,pt, & 1481 reffrad,profmmr) 1482 zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_n2)) & ! AF: TODO: replace by igcm_haze? 1483 /ptimestep 1484 else 1485 !! AF: TODO import from pluto.old? 1486 ! call hazecloud(ngrid,nlayer,nq,ptimestep,& 1487 ! pplay,pplev,pq,pdq,dist_sol,mu0,zfluxuv,zdqhaze,& 1488 ! zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin) 1489 endif 1490 1491 DO iq=1, nq ! should be updated 1492 DO l=1,nlayer 1493 DO ig=1,ngrid 1494 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqhaze(ig,l,iq) 1495 ENDDO 1496 ENDDO 1497 ENDDO 1422 1498 endif ! end of 'tracer' 1423 1499 … … 1431 1507 1432 1508 call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic) 1433 1509 1434 1510 if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer 1435 1511 1436 1512 do l = 1, nlayer 1437 1513 do ig=1,ngrid … … 1442 1518 1443 1519 end if 1444 1520 1445 1521 end do ! iq=1,nq 1446 1522 … … 1449 1525 1450 1526 if (is_master) print*,'--> Ls =',zls*180./pi 1451 1452 1527 1528 1453 1529 !---------------------------------------------------------------------- 1454 1530 ! Writing NetCDF file "RESTARTFI" at the end of the run … … 1468 1544 !! Update surface ice distribution to iterate to steady state if requested 1469 1545 !! AF24: removed 1470 1546 1471 1547 ! endif 1472 1548 #ifndef MESOSCALE 1473 1549 1474 1550 if (ngrid.ne.1) then 1475 1551 write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin 1476 1552 1477 1553 call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, & 1478 1554 ptimestep,ztime_fin, & … … 1493 1569 ! which can later be used to make the statistic files of the run: 1494 1570 ! if flag "callstats" from callphys.def is .true.) 1495 1571 1496 1572 1497 1573 call wstats(ngrid,"ps","Surface pressure","Pa",2,ps) … … 1503 1579 "Thermal IR radiative flux to space","W.m-2",2, & 1504 1580 fluxtop_lw) 1505 1581 1506 1582 ! call wstats(ngrid,"fluxsurf_sw", & 1507 1583 ! "Solar radiative flux to surface","W.m-2",2, & 1508 ! fluxsurf_sw_tot) 1584 ! fluxsurf_sw_tot) 1509 1585 ! call wstats(ngrid,"fluxtop_sw", & 1510 1586 ! "Solar radiative flux to space","W.m-2",2, & … … 1527 1603 do iq=1,nq 1528 1604 call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq)) 1529 call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & 1605 call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & 1530 1606 'kg m^-2',2,qsurf(1,iq) ) 1531 call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & 1607 call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & 1532 1608 'kg m^-2',2,qcol(1,iq) ) 1533 1534 ! call wstats(ngrid,trim(noms(iq))//'_reff', & 1535 ! trim(noms(iq))//'_reff', & 1609 1610 ! call wstats(ngrid,trim(noms(iq))//'_reff', & 1611 ! trim(noms(iq))//'_reff', & 1536 1612 ! 'm',3,reffrad(1,1,iq)) 1537 1613 1538 1614 end do 1539 1540 endif ! end of 'tracer' 1615 1616 endif ! end of 'tracer' 1541 1617 1542 1618 !AF24: deleted slab ocean and water … … 1546 1622 call mkstats(ierr) 1547 1623 endif 1548 1624 1549 1625 1550 1626 #ifndef MESOSCALE 1551 1627 1552 1628 !----------------------------------------------------------------------------------------------------- 1553 1629 ! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic … … 1577 1653 1578 1654 ! ! Oceanic layers !AF24: removed 1579 1655 1580 1656 ! ! Thermal plume model !AF24: removed 1581 1657 1582 1658 ! GW non-oro outputs !AF24: removed 1583 1659 1584 1660 ! Total energy balance diagnostics 1585 1661 if(callrad)then 1586 1662 1587 1663 !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent) 1588 1664 !call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1)) … … 1591 1667 call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw) 1592 1668 call writediagfi(ngrid,"shad","rings"," ", 2, fract) 1593 1669 1594 1670 ! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1) 1595 1671 ! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1) … … 1604 1680 call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd) 1605 1681 ! endif 1606 1682 1607 1683 call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn) 1608 1684 1609 1685 endif ! end of 'callrad' 1610 1686 1611 1687 if(enertest) then 1612 1688 1613 1689 if (calldifv) then 1614 1690 1615 1691 call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2) 1616 1692 call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux) 1617 1693 1618 1694 ! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff) 1619 1695 ! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff) 1620 1696 ! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs) 1621 1622 endif 1623 1697 1698 endif 1699 1624 1700 if (corrk) then 1625 1701 call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw) 1626 1702 call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw) 1627 1703 endif 1628 1704 1629 1705 ! if(watercond) then !AF24: removed 1630 1706 … … 1633 1709 call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE) 1634 1710 call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation) 1635 1636 endif 1637 1711 1712 endif 1713 1638 1714 endif ! end of 'enertest' 1639 1715 1640 1716 ! Diagnostics of optical thickness 1641 1717 ! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19 1642 if (diagdtau) then 1718 if (diagdtau) then 1643 1719 do nw=1,L_NSPECTV 1644 1720 write(str2,'(i2.2)') nw … … 1657 1733 call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad) 1658 1734 call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn) 1659 1735 1660 1736 ! For Debugging. 1661 1737 !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat)) 1662 1738 !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi) 1663 1664 1665 ! Output aerosols. 1666 if (igcm_n2_ice.ne.0.and.iaero_n2.ne.0) &1667 call writediagfi(ngrid,'N2ice_reff','N2ice_reff','m',3,reffrad(1,1,iaero_n2))1668 if (igcm_n2_ice.ne.0.and.iaero_n2.ne.0) &1669 call writediagfi(ngrid,'N2ice_reffcol','N2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_n2))1739 1740 1741 ! Output aerosols.!AF: TODO: write haze aerosols 1742 ! if (igcm_n2_ice.ne.0.and.iaero_haze.ne.0) & 1743 ! call writediagfi(ngrid,'N2ice_reff','N2ice_reff','m',3,reffrad(1,1,iaero_haze)) 1744 ! if (igcm_n2_ice.ne.0.and.iaero_haze.ne.0) & 1745 ! call writediagfi(ngrid,'N2ice_reffcol','N2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_haze)) 1670 1746 ! if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) & !AF24: removed 1671 1747 1672 1748 ! Output tracers. 1673 1749 if (tracer) then 1674 1750 1675 1751 do iq=1,nq 1676 1752 call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq)) 1677 call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & 1753 call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & 1678 1754 'kg m^-2',2,qsurf_hist(1,iq) ) 1679 call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & 1755 call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & 1680 1756 'kg m^-2',2,qcol(1,iq) ) 1681 ! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & 1682 ! 'kg m^-2',2,qsurf(1,iq) ) 1757 ! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & 1758 ! 'kg m^-2',2,qsurf(1,iq) ) 1683 1759 1684 1760 ! if(watercond.or.CLFvarying)then !AF24: removed 1685 1761 1686 1762 if(generic_condensation)then 1687 1763 call writediagfi(ngrid,"rneb_generic","GCS cloud fraction (generic condensation)"," ",3,rneb_generic) … … 1692 1768 ! if(generic_rain)then !AF24: removed 1693 1769 ! if((hydrology).and.(.not.ok_slab_ocean))then !AF24: removed 1694 1770 1695 1771 call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col) 1696 1772 1697 1773 enddo ! end of 'nq' loop 1698 1774 1699 1775 endif ! end of 'tracer' 1700 1776 1701 1777 1702 1778 ! Output spectrum. 1703 if(specOLR.and.corrk)then 1779 if(specOLR.and.corrk)then 1704 1780 call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu) 1705 1781 call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu) … … 1713 1789 if (.not.calldifv) comm_LATENT_HF(:)=0.0 1714 1790 ! if ((tracer).and.(water)) then !AF24: removed 1715 1791 1716 1792 if ((tracer).and.(generic_condensation)) then 1717 1793 ! .and.(.not. water) … … 1724 1800 do iq=1,nq 1725 1801 call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic) 1726 1802 1727 1803 if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer 1728 1804 … … 1741 1817 comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq) 1742 1818 1743 endif 1744 end do ! do iq=1,nq loop on tracers 1819 endif 1820 end do ! do iq=1,nq loop on tracers 1745 1821 1746 1822 else … … 1785 1861 1786 1862 ! XIOS outputs 1787 #ifdef CPP_XIOS 1863 #ifdef CPP_XIOS 1788 1864 ! Send fields to XIOS: (NB these fields must also be defined as 1789 1865 ! <field id="..." /> in context_lmdz_physics.xml to be correctly used) 1790 1866 CALL send_xios_field("ls",zls) 1791 1867 1792 1868 CALL send_xios_field("ps",ps) 1793 1869 CALL send_xios_field("area",cell_area) … … 1797 1873 CALL send_xios_field("v",zv) 1798 1874 CALL send_xios_field("omega",omega) 1799 1875 1800 1876 ! IF (calltherm) THEN !AF24: removed 1801 1877 ! IF (water) THEN !AF24: removed 1802 1878 1803 1879 CALL send_xios_field("ISR",fluxtop_dn) 1804 1880 CALL send_xios_field("OLR",fluxtop_lw) 1805 1881 CALL send_xios_field("ASR",fluxabs_sw) 1806 1882 1807 if (specOLR .and. corrk) then 1883 if (specOLR .and. corrk) then 1808 1884 call send_xios_field("OLR3D",OLR_nu) 1809 1885 call send_xios_field("IR_Bandwidth",DWNI) … … 1818 1894 endif 1819 1895 #endif 1820 1896 1821 1897 if (check_physics_outputs) then 1822 1898 ! Check the validity of updated fields at the end of the physics step … … 1825 1901 1826 1902 icount=icount+1 1827 1903 1828 1904 end subroutine physiq 1829 1905 1830 1906 end module physiq_mod -
trunk/LMDZ.PLUTO/libf/phypluto/radii_mod.F90
r3184 r3195 4 4 ! module to centralize the radii calculations for aerosols 5 5 !================================================================== 6 6 7 7 ! N2 cloud properties (initialized in inifis) 8 8 real,save :: Nmix_n2 ! Number mixing ratio of N2 ice particles … … 21 21 ! Purpose 22 22 ! ------- 23 ! Compute the effective radii of liquid and icy water particles 24 ! Jeremy Leconte (2012) 25 ! Extended to dust, N2, NH3, 2-lay,Nlay,auroral aerosols by ?? 26 ! Added Radiative Generic Condensable Species effective radii 27 ! calculations (Lucas Teinturier 2022) 23 ! Compute the effective radii of haze particles (TB) 28 24 ! 29 ! Authors30 ! -------31 ! Jeremy Leconte (2012)32 25 ! 33 26 !================================================================== … … 35 28 use ioipsl_getin_p_mod, only: getin_p 36 29 use radinc_h, only: naerkind 37 use aerosol_mod, only: iaero_back2lay, iaero_n2, iaero_dust, & 38 iaero_h2so4, iaero_nh3, iaero_nlay, & 39 iaero_aurora, iaero_generic, i_rgcs_ice 40 use callkeys_mod, only: size_nh3_cloud, nlayaero, aeronlay_size, & 30 use aerosol_mod, only: iaero_haze, i_haze, & 31 iaero_generic, i_rgcs_ice 32 use callkeys_mod, only: nlayaero, aeronlay_size, & 41 33 aeronlay_nueff,aerogeneric 42 use tracer_h, only: radius, nqtot, is_rgcs 34 use tracer_h, only: radius, nqtot, is_rgcs, nmono 43 35 Implicit none 44 36 … … 47 39 48 40 real, intent(out) :: reffrad(ngrid,nlayer,naerkind) !aerosols radii (K) 49 real, intent(out) :: nueffrad(ngrid,nlayer,naerkind) !variance 41 real, intent(out) :: nueffrad(ngrid,nlayer,naerkind) !variance 50 42 51 logical, save :: firstcall=.true. 52 !$OMP THREADPRIVATE(firstcall) 53 integer :: iaer, ia , iq, i_rad 43 integer :: iaer, ia , iq, i_rad 54 44 55 45 do iaer=1,naerkind 56 ! these values will change once the microphysics gets to work 57 ! UNLESS tracer=.false., in which case we should be working with 58 ! a fixed aerosol layer, and be able to define reffrad in a 59 ! .def file. To be improved! 60 ! |-> Done in th n-layer aerosol case (JVO 20) 61 62 if(iaer.eq.iaero_n2)then ! N2 ice 63 reffrad(1:ngrid,1:nlayer,iaer) = 1.e-4 64 nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 46 print*, 'TB22 iaer',iaer,iaero_haze 47 if(iaer.eq.iaero_haze)then 48 ! Equivalent sphere radius 49 reffrad(1:ngrid,1:nlayer,iaer)=radius(i_haze)*nmono**(1./3.) 50 !reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6 ! haze 51 nueffrad(1:ngrid,1:nlayer,iaer) = 0.02 ! haze 52 print*, 'TB22 Hello2',radius(i_haze)*nmono**(1./3.) 65 53 endif 66 67 if(iaer.eq.iaero_dust)then ! dust 68 reffrad(1:ngrid,1:nlayer,iaer) = 1.e-5 69 nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 70 endif 71 72 if(iaer.eq.iaero_h2so4)then ! H2SO4 ice 73 reffrad(1:ngrid,1:nlayer,iaer) = 1.e-6 74 nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 75 endif 76 77 if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols 78 reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6 79 nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 80 endif 81 82 83 if(iaer.eq.iaero_nh3)then ! Nh3 cloud 84 reffrad(1:ngrid,1:nlayer,iaer) = size_nh3_cloud 85 nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 86 endif 87 88 do ia=1,nlayaero 89 if(iaer.eq.iaero_nlay(ia))then ! N-layer aerosols 90 reffrad(1:ngrid,1:nlayer,iaer) = aeronlay_size(ia) 91 nueffrad(1:ngrid,1:nlayer,iaer) = aeronlay_nueff(ia) 92 endif 93 enddo 94 95 if(iaer.eq.iaero_aurora)then ! Auroral aerosols 96 reffrad(1:ngrid,1:nlayer,iaer) = 3.e-7 97 nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 98 endif 99 100 do ia=1,aerogeneric ! Radiative Generic Condensable Species 101 if (iaer .eq. iaero_generic(ia)) then 102 i_rad = i_rgcs_ice(ia) 103 reffrad(1:ngrid,1:nlayer,iaer)=radius(i_rad) 104 nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 105 endif 106 enddo ! generic radiative condensable aerosols 107 108 enddo ! iaer=1,naerkind 54 enddo 109 55 110 56 end subroutine su_aer_radii … … 113 59 114 60 !================================================================== 115 subroutine n2_reffrad(ngrid,nlayer,nq,pq,reffrad)116 !==================================================================61 subroutine haze_reffrad_fix(ngrid,nlayer,zzlay,reffrad,nueffrad) 62 !================================================================== 117 63 ! Purpose 118 64 ! ------- 119 ! Compute the effective radii of n2 ice particles !AF24: copied from CO2 120 ! 121 ! Authors 122 ! ------- 123 ! Jeremy Leconte (2012) 65 ! Compute the effective radii of haze particles 66 ! fixed profile of radius (TB) 124 67 ! 125 68 !================================================================== 126 USE tracer_h, only:igcm_n2_ice,rho_n2 69 use radinc_h, only: naerkind 70 USE tracer_h, only:rho_n2,nmono 127 71 use comcstfi_mod, only: pi 72 use aerosol_mod, only: iaero_haze, i_haze 73 use datafile_mod 128 74 Implicit none 129 75 130 integer,intent(in) :: ngrid,nlayer,nq 76 integer,intent(in) :: ngrid,nlayer 77 real,intent(in) :: zzlay(ngrid,nlayer) 78 real, intent(out) :: reffrad(ngrid,nlayer,naerkind) ! haze particles radii (m) 79 real, intent(out) :: nueffrad(ngrid,nlayer,naerkind) ! 131 80 132 real , intent(in) :: pq(ngrid,nlayer,nq) !tracer mixing ratios (kg/kg)133 real, intent(out) :: reffrad(ngrid,nlayer) !n2 ice particles radii (m)81 real :: zrad 82 real,external :: CBRT 134 83 135 integer :: ig,l 136 real :: zrad 137 real,external :: CBRT 138 139 84 logical, save :: firstcall=.true. 85 !$OMP THREADPRIVATE(firstcall) 140 86 141 if (radfixed) then 142 reffrad(1:ngrid,1:nlayer) = 5.e-5 ! N2 ice 143 else 144 do l=1,nlayer 145 do ig=1,ngrid 146 zrad = CBRT( 3*pq(ig,l,igcm_n2_ice)/(4*Nmix_n2*pi*rho_n2) ) 147 reffrad(ig,l) = min(max(zrad,1.e-6),100.e-6) 148 enddo 149 enddo 150 end if 87 ! Local variables 88 integer :: iaer,l,ifine,ig 89 real :: radvec(ngrid,nlayer) 151 90 152 end subroutine n2_reffrad 91 !!read altitudes and radius 92 integer Nfine 93 !parameter(Nfine=21) 94 parameter(Nfine=701) 95 character(len=100) :: file_path 96 real,save :: levdat(Nfine),raddat(Nfine) 97 98 !---------------- INPUT ------------------------------------------------ 99 100 IF (firstcall) then 101 firstcall=.false. 102 file_path=trim(datadir)//'/haze_prop/hazerad.txt' 103 open(223,file=file_path,form='formatted') 104 do ifine=1,Nfine 105 read(223,*) levdat(ifine), raddat(ifine) 106 enddo 107 close(223) 108 print*, 'TB22 READ HAZERAD' 109 ENDIF 110 111 ! in radii mod levs has been put in km 112 DO ig=1,ngrid 113 CALL interp_line(levdat,raddat,Nfine,zzlay(ig,:)/1000,radvec(ig,:),nlayer) 114 enddo 115 116 do iaer=1,naerkind 117 if(iaer.eq.iaero_haze)then 118 ! spherical radius or eq spherical radius 119 ! TB22: fractal has no impact on reffrad if haze_radproffix 120 do ig=1,ngrid 121 do l=1,nlayer 122 reffrad(ig,l,iaer)=radvec(ig,l)*1.e-9 ! nm => m 123 enddo 124 enddo 125 nueffrad(1:ngrid,1:nlayer,iaer) = 0.02 ! haze 126 endif 127 enddo 128 129 end subroutine haze_reffrad_fix 153 130 !================================================================== 154 131 155 132 156 157 !==================================================================158 subroutine dust_reffrad(ngrid,nlayer,reffrad)159 !==================================================================160 ! Purpose161 ! -------162 ! Compute the effective radii of dust particles163 !164 ! Authors165 ! -------166 ! Jeremy Leconte (2012)167 !168 !==================================================================169 Implicit none170 171 integer,intent(in) :: ngrid172 integer,intent(in) :: nlayer173 174 real, intent(out) :: reffrad(ngrid,nlayer) !dust particles radii (m)175 176 reffrad(1:ngrid,1:nlayer) = 2.e-6 ! dust177 178 end subroutine dust_reffrad179 !==================================================================180 181 182 !==================================================================183 subroutine h2so4_reffrad(ngrid,nlayer,reffrad)184 !==================================================================185 ! Purpose186 ! -------187 ! Compute the effective radii of h2so4 particles188 !189 ! Authors190 ! -------191 ! Jeremy Leconte (2012)192 !193 !==================================================================194 Implicit none195 196 integer,intent(in) :: ngrid197 integer,intent(in) :: nlayer198 199 real, intent(out) :: reffrad(ngrid,nlayer) !h2so4 particle radii (m)200 201 reffrad(1:ngrid,1:nlayer) = 1.e-6 ! h2so4202 203 end subroutine h2so4_reffrad204 !==================================================================205 206 !==================================================================207 subroutine back2lay_reffrad(ngrid,reffrad,nlayer,pplev)208 !==================================================================209 ! Purpose210 ! -------211 ! Compute the effective radii of particles in a 2-layer model212 !213 ! Authors214 ! -------215 ! Sandrine Guerlet (2013)216 !217 !==================================================================218 use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo,size_tropo, &219 pres_bottom_strato,size_strato220 221 Implicit none222 223 integer,intent(in) :: ngrid224 225 real, intent(out) :: reffrad(ngrid,nlayer) ! particle radii (m)226 REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)227 INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers228 REAL :: expfactor229 INTEGER l,ig230 231 reffrad(:,:)=1e-6 !!initialization, not important232 DO ig=1,ngrid233 DO l=1,nlayer-1234 IF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN235 reffrad(ig,l) = size_tropo236 ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN237 expfactor=log(size_strato/size_tropo) / log(pres_bottom_strato/pres_top_tropo)238 reffrad(ig,l)= size_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)239 ELSEIF (pplev(ig,l) .le. pres_bottom_strato) then240 reffrad(ig,l) = size_strato241 ENDIF242 ENDDO243 ENDDO244 245 end subroutine back2lay_reffrad246 !==================================================================247 248 133 end module radii_mod 249 134 !================================================================== -
trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90
r3193 r3195 16 16 use radcommon_h, only: radiustab,nsize,tstellar 17 17 use radcommon_h, only: qrefvis,qrefir,omegarefir !,omegarefvis 18 use aerosol_mod, only: noaero,iaero_n2,iaero_dust,iaero_h2so4 19 use aerosol_mod, only: iaero_back2lay,iaero_nh3,iaero_nlay,iaero_aurora 20 use aerosol_mod, only: iaero_generic,i_rgcs_ice 18 use aerosol_mod, only: noaero,iaero_haze 21 19 use callkeys_mod, only: tplanet, optprop_back2lay_vis, optprop_back2lay_ir, & 22 20 optprop_aeronlay_vis, optprop_aeronlay_ir, & 23 21 aeronlay_lamref, nlayaero,aerogeneric 24 22 use tracer_h, only: noms 25 23 26 24 use mod_phys_lmdz_para, only : is_master, bcast 27 25 … … 41 39 ! omeg(nwvl) omegav omegavis(L_NSPECTV) 42 40 ! gfactor(nwvl) gav gvis(L_NSPECTV) 43 ! 41 ! 44 42 ! Authors 45 ! ------- 43 ! ------- 46 44 ! Richard Fournier (1996) Francois Forget (1996) 47 45 ! Frederic Hourdin … … 54 52 ! - Add water ice clouds 55 53 ! MODIF R. Wordsworth (2009) 56 ! - generalisation to arbitrary spectral bands 54 ! - generalisation to arbitrary spectral bands 57 55 ! 58 56 !================================================================== … … 100 98 REAL omegavIR(L_NSPECTI) ! Average single scattering albedo 101 99 REAL gavIR(L_NSPECTI) ! Average assymetry parameter 102 100 103 101 logical forgetit ! use francois' data? 104 102 integer iwvl, ia … … 109 107 110 108 !---- Please indicate the names of the optical property files below 111 ! Please also choose the reference wavelengths of each aerosol 109 ! Please also choose the reference wavelengths of each aerosol 112 110 113 111 !--------- README TO UNDERSTAND WHAT FOLLOWS (JVO 20) ------- … … 136 134 print*, 'naerkind= 0' 137 135 endif 138 139 136 137 140 138 do iaer=1,naerkind 141 139 ! naerkind=1: haze 142 140 if (iaer.eq.1) then 143 141 144 ! Only one table of optical properties : 142 ! Only one table of optical properties : 145 143 write(*,*)'Suaer haze optical properties, using: ', & 146 144 TRIM(hazeprop) … … 150 148 ! infrared 151 149 file_id(iaer,2)=file_id(iaer,1) 152 150 153 151 lamrefvis(iaer)=0.607E-6 ! reference wavelength for opacity vis pivot wavelength Cheng et al 2017 154 152 lamrefir(iaer)=2.E-6 ! reference wavelength for opacity IR (in the LEISA range) 155 153 156 154 endif 157 enddo 155 enddo 158 156 159 157 IF (naerkind .gt. 1) THEN … … 162 160 call abort 163 161 ENDIF 164 162 165 163 !------------------------------------------------------------------ 166 164 … … 185 183 !!!!$OMP MASTER 186 184 if (is_master) then 187 185 188 186 INQUIRE(FILE=TRIM(datadir)//'/'//TRIM(aerdir)//& 189 187 '/'//TRIM(file_id(iaer,idomain)),& … … 205 203 '/'//TRIM(file_id(iaer,idomain)),& 206 204 FORM='formatted') 207 ENDIF 205 ENDIF 208 206 ENDIF 209 207 IF(.NOT.file_ok) THEN … … 220 218 write(*,*)' http://www.lmd.jussieu.fr/',& 221 219 '~lmdz/planets/LMDZ.GENERIC/datagcm/' 222 CALL ABORT 220 CALL ABORT 223 221 ENDIF 224 222 … … 239 237 read(file_unit,*) nsize(iaer,idomain) 240 238 endwhile = .true. 241 CASE DEFAULT reading1_seq ! ---------------------------- 242 WRITE(*,*) 'readoptprop: ',& 243 'Error while loading optical properties.' 244 CALL ABORT 239 CASE DEFAULT reading1_seq ! ---------------------------- 240 WRITE(*,*) 'readoptprop: ',& 241 'Error while loading optical properties.' 242 CALL ABORT 245 243 END SELECT reading1_seq ! ============================== 246 244 ENDIF … … 256 254 ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn 257 255 ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep 258 ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg 256 ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg 259 257 ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g 260 258 … … 316 314 317 315 jfile = jfile+1 318 IF ((idomain.NE.iaero_ n2).OR.(iaer.NE.iaero_n2)) THEN316 IF ((idomain.NE.iaero_haze).OR.(iaer.NE.iaero_haze)) THEN 319 317 endwhile = .true. 320 318 ENDIF … … 334 332 335 333 ! 1.5 If Francois' data, convert wvl to metres 336 if(iaer.eq.iaero_ n2)then334 if(iaer.eq.iaero_haze)then 337 335 if(forgetit)then 338 336 ! print*,'please sort out forgetit for naerkind>1' … … 440 438 DEALLOCATE(wvl) ! wvl 441 439 DEALLOCATE(radiusdyn) ! radiusdyn 442 DEALLOCATE(ep) ! ep 443 DEALLOCATE(omeg) ! omeg 440 DEALLOCATE(ep) ! ep 441 DEALLOCATE(omeg) ! omeg 444 442 DEALLOCATE(gfactor) ! g 445 443 … … 447 445 END DO ! Loop on idomain 448 446 !======================================================================== 449 447 450 448 ! cleanup 451 449 deallocate(file_id) 452 450 453 451 end subroutine suaer_corrk 454 452 455 453 end module suaer_corrk_mod -
trunk/LMDZ.PLUTO/libf/phypluto/surfdat_h.F90
r3184 r3195 16 16 real,allocatable,dimension(:) :: zmea,zstd,zsig,zgam,zthe 17 17 !$OMP THREADPRIVATE(zmea,zstd,zsig,zgam,zthe) 18 real ttop 19 real,allocatable,dimension(:) :: kp ! TB ref pressure 20 real p00 21 !$OMP THREADPRIVATE(ttop,kp,p00) 22 ! surface properties ! TB16 23 real alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a 24 !$OMP THREADPRIVATE(alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a) 25 real emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq 26 !$OMP THREADPRIVATE(emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq) 27 real ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe 28 !$OMP THREADPRIVATE(ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe) 29 real alb_tho_spe,emis_tho_spe 30 !$OMP THREADPRIVATE(alb_tho_spe,emis_tho_spe) 18 31 19 32 real,allocatable,dimension(:) :: dryness !"Dryness coefficient" for grnd water ice sublimation -
trunk/LMDZ.PLUTO/libf/phypluto/surfini.F
r3184 r3195 3 3 4 4 USE surfdat_h, only: albedodat 5 USE tracer_h, only: igcm_ n2_ice5 USE tracer_h, only: igcm_haze 6 6 use planetwide_mod, only: planetwide_maxval, planetwide_minval 7 7 use radinc_h, only : L_NSPECTV … … 9 9 10 10 IMPLICIT NONE 11 12 11 12 13 13 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 14 14 cccccccccccccc cccccccccccccc 15 15 cccccccccccccc Spectral Albedo Initialisation - Routine modified by MT2015. cccccccccccccc 16 cccccccccccccc cccccccccccccc 16 cccccccccccccc cccccccccccccc 17 17 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 18 18 … … 56 56 57 57 ! Step 3 : We modify the albedo considering some N2 at the surface. We dont take into account water ice (this is made in hydrol after the first timestep) ... 58 if (igcm_ n2_ice.ne.0) then58 if (igcm_haze.ne.0) then 59 59 DO ig=1,ngrid 60 IF (qsurf(ig,igcm_ n2_ice) .GT. 1.) THEN ! This was changed by MT2015. Condition for ~1mm of N2 ice deposit.60 IF (qsurf(ig,igcm_haze) .GT. 1.) THEN ! This was changed by MT2015. Condition for ~1mm of N2 ice deposit. 61 61 DO nw=1,L_NSPECTV 62 62 albedo(ig,nw)=albedo_n2_ice_SPECTV(nw) 63 63 ENDDO 64 END IF 65 ENDDO 64 END IF 65 ENDDO 66 66 else 67 67 write(*,*) "surfini: No N2 ice tracer on surface ..." 68 68 write(*,*) " and therefore no albedo change." 69 endif 69 endif 70 70 call planetwide_minval(albedo,min_albedo) 71 71 call planetwide_maxval(albedo,max_albedo) -
trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90
r3184 r3195 14 14 15 15 character*30, save, allocatable :: noms(:) ! name of the tracer 16 real, save, allocatable :: mmol(:) ! mole mass of tracer (g/mol) 16 real, save, allocatable :: mmol(:) ! mole mass of tracer (g/mol) 17 17 real, save, allocatable :: aki(:) ! to compute coefficient of thermal concduction if photochem 18 18 real, save, allocatable :: cpi(:) ! to compute cpnew in concentration.F if photochem … … 28 28 real,save :: rho_dust ! Mars dust density (kg.m-3) 29 29 real,save :: rho_ice ! Water ice density (kg.m-3) 30 real,save :: rho_ch4_ice ! ch4 ice density (kg.m-3) 31 real,save :: rho_co_ice ! co ice density (kg.m-3) 30 32 real,save :: rho_n2 ! N2 ice density (kg.m-3) 33 real,save :: lw_ch4 ! Latent heat CH4 gas -> solid 34 real,save :: lw_co ! Latent heat CO gas -> solid 35 real,save :: lw_n2 ! Latent heat N2 gas -> solid 36 integer,save :: nmono 31 37 real,save :: ref_r0 ! for computing reff=ref_r0*r0 (in log.n. distribution) 32 38 !$OMP THREADPRIVATE(noms,mmol,aki,cpi,radius,rho_q,qext,alpha_lift,alpha_devil,qextrhor, & 33 !$OMP varian,r3n_q,rho_dust,rho_ice,rho_n2, ref_r0)39 !$OMP varian,r3n_q,rho_dust,rho_ice,rho_n2,lw_n2,ref_r0) 34 40 35 41 integer, save, allocatable :: is_chim(:) ! 1 if tracer used in chemistry, else 0 … … 59 65 ! tracer indexes: these are initialized in initracer and should be 0 if the 60 66 ! corresponding tracer does not exist 61 ! dust 62 integer,save,allocatable :: igcm_dustbin(:) ! for dustbin 'dust' tracers 63 ! dust, special doubleq case 64 integer,save :: igcm_dust_mass ! dust mass mixing ratio (for transported dust) 65 integer,save :: igcm_dust_number ! dust number mixing ratio (transported dust) 66 ! water 67 integer,save :: igcm_h2o_vap ! water vapour 68 integer,save :: igcm_h2o_ice ! water ice 69 ! chemistry: 70 integer,save :: igcm_co2 71 integer,save :: igcm_co 72 integer,save :: igcm_o 73 integer,save :: igcm_o1d 74 integer,save :: igcm_o2 75 integer,save :: igcm_o3 76 integer,save :: igcm_h 77 integer,save :: igcm_h2 78 integer,save :: igcm_oh 79 integer,save :: igcm_ho2 80 integer,save :: igcm_h2o2 67 68 !Pluto chemistry 69 integer,save :: igcm_co_gas 81 70 integer,save :: igcm_n2 82 71 integer,save :: igcm_ar 83 integer,save :: igcm_n 84 integer,save :: igcm_no 85 integer,save :: igcm_no2 86 integer,save :: igcm_n2d 87 integer,save :: igcm_ch4 72 integer,save :: igcm_ch4_gas ! methane gas 73 ! other tracers 74 integer,save :: igcm_ar_n2 ! for simulations using co2 +neutral gaz 75 integer,save :: igcm_ch4_ice ! methane ice 76 integer,save :: igcm_co_ice ! methane ice 77 integer,save :: igcm_prec_haze 78 integer,save :: igcm_haze 79 integer,save :: igcm_haze10 80 integer,save :: igcm_haze30 81 integer,save :: igcm_haze50 82 integer,save :: igcm_haze100 83 integer,save :: igcm_eddy1e6 84 integer,save :: igcm_eddy1e7 85 integer,save :: igcm_eddy5e7 86 integer,save :: igcm_eddy1e8 87 integer,save :: igcm_eddy5e8 88 88 89 integer,save :: igcm_ch3 90 integer,save :: igcm_ch 91 integer,save :: igcm_3ch2 92 integer,save :: igcm_1ch2 93 integer,save :: igcm_cho 94 integer,save :: igcm_ch2o 95 integer,save :: igcm_ch3o 96 integer,save :: igcm_c 97 integer,save :: igcm_c2 98 integer,save :: igcm_c2h 99 integer,save :: igcm_c2h2 100 integer,save :: igcm_c2h3 101 integer,save :: igcm_c2h4 102 integer,save :: igcm_c2h6 103 integer,save :: igcm_ch2co 104 integer,save :: igcm_ch3co 105 integer,save :: igcm_hcaer 106 107 ! other tracers 108 integer,save :: igcm_ar_n2 ! for simulations using n2 +neutral gaz 109 integer,save :: igcm_n2_ice ! N2 ice 110 !$OMP THREADPRIVATE(igcm_dustbin,igcm_dust_mass,igcm_dust_number,igcm_h2o_vap,igcm_h2o_ice, & 111 !$OMP igcm_co2,igcm_co,igcm_o,igcm_o1d,igcm_o2,igcm_o3,igcm_h,igcm_h2,igcm_oh, & 112 !$OMP igcm_ho2,igcm_h2o2,igcm_n2,igcm_ar,igcm_ar_n2,igcm_n2_ice, & 113 !$OMP igcm_n,igcm_no,igcm_no2,igcm_n2d,igcm_ch4,igcm_ch3,igcm_ch,igcm_3ch2, & 114 !$OMP igcm_1ch2,igcm_cho,igcm_ch2o,igcm_ch3o,igcm_c,igcm_c2,igcm_c2h,igcm_c2h2, & 115 !$OMP igcm_c2h3,igcm_c2h4,igcm_c2h6,igcm_ch2co,igcm_ch3co,igcm_hcaer) 89 !$OMP THREADPRIVATE(igcm_co_gas,igcm_n2,igcm_ar,igcm_ch4_gas,igcm_ar_n2,&igcm_ch4_ice,igcm_co_ice,igcm_prec_haze,igcm_haze,igcm_haze10,igcm_haze30,igcm_haze50,igcm_haze100,igcm_eddy1e6,igcm_eddy1e7,igcm_eddy5e7,igcm_eddy1e8,igcm_eddy5e8) 116 90 117 91 end module tracer_h -
trunk/LMDZ.PLUTO/libf/phypluto/vdifc_mod.F
r3184 r3195 1 1 module vdifc_mod 2 2 3 3 implicit none 4 4 5 5 contains 6 7 subroutine vdifc(ngrid,nlay,nq,ppopsk, 8 & ptimestep,pcapcal,lecrit, 6 7 subroutine vdifc(ngrid,nlay,nq,ppopsk, 8 & ptimestep,pcapcal,lecrit, 9 9 & pplay,pplev,pzlay,pzlev,pz0, 10 10 & pu,pv,ph,pq,ptsrf,pemis,pqsurf, … … 13 13 & pdqdif,pdqsdif) 14 14 15 ! use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol 15 ! use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol 16 16 ! & ,Psat_water, Lcpdqsat_water 17 17 use radcommon_h, only : sigma 18 18 use surfdat_h, only: dryness 19 use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice20 19 use comcstfi_mod, only: g, r, cpp, rcp 21 20 use callkeys_mod, only: tracer,nosurf … … 24 23 25 24 !================================================================== 26 ! 25 ! 27 26 ! Purpose 28 27 ! ------- 29 28 ! Turbulent diffusion (mixing) for pot. T, U, V and tracers 30 ! 29 ! 31 30 ! Implicit scheme 32 31 ! We start by adding to variables x the physical tendencies … … 34 33 ! 35 34 ! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1) 36 ! 35 ! 37 36 ! Authors 38 ! ------- 37 ! ------- 39 38 ! F. Hourdin, F. Forget, R. Fournier (199X) 40 39 ! R. Wordsworth, B. Charnay (2010) 41 ! 40 ! 42 41 !================================================================== 43 42 … … 62 61 REAL,INTENT(IN) :: pcapcal(ngrid) 63 62 REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1) 64 63 65 64 ! Arguments added for condensation 66 65 REAL,INTENT(IN) :: ppopsk(ngrid,nlay) … … 70 69 ! Tracers 71 70 ! -------- 72 integer,intent(in) :: nq 71 integer,intent(in) :: nq 73 72 real,intent(in) :: pqsurf(ngrid,nq) 74 real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 75 real,intent(out) :: pdqdif(ngrid,nlay,nq) 76 real,intent(out) :: pdqsdif(ngrid,nq) 77 73 real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 74 real,intent(out) :: pdqdif(ngrid,nlay,nq) 75 real,intent(out) :: pdqsdif(ngrid,nq) 76 78 77 ! local 79 78 ! ----- … … 101 100 SAVE firstcall 102 101 !$OMP THREADPRIVATE(firstcall) 103 102 104 103 ! variables added for N2 condensation 105 104 ! ------------------------------------ … … 134 133 real masse, Wtot, Wdiff 135 134 136 real dqsdif_total(ngrid) 137 real zq0(ngrid) 135 real dqsdif_total(ngrid) 136 real zq0(ngrid) 138 137 139 138 forceWC=.true. … … 151 150 ! PRINT*,' acond,bcond',acond,bcond 152 151 153 ! if(water)then154 ! iliq=igcm_h2o_vap155 ! ivap=igcm_h2o_vap156 ! iice=igcm_h2o_ice ! simply to make the code legible157 ! ! to be generalised later158 ! endif159 152 160 153 firstcall=.false. … … 208 201 DO ilev=1,nlay 209 202 DO ig=1,ngrid 210 zq(ig,ilev,iq)=pq(ig,ilev,iq) + 203 zq(ig,ilev,iq)=pq(ig,ilev,iq) + 211 204 & pdqfi(ig,ilev,iq)*ptimestep 212 205 ENDDO … … 220 213 ! 221 214 ! Source of turbulent kinetic energy at the surface 222 ! ------------------------------------------------- 215 ! ------------------------------------------------- 223 216 ! Formula is Cd_0 = (karman / log[1+z1/z0])^2 224 217 … … 242 235 243 236 ! Turbulent diffusion coefficients in the boundary layer 244 ! ------------------------------------------------------ 237 ! ------------------------------------------------------ 245 238 246 239 call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay … … 272 265 ! donc les entrees sont /zcdv/ pour la condition a la limite sol 273 266 ! et /zkv/ = Ku 274 267 275 268 CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2)) 276 269 CALL multipl(ngrid,zcdv,zb0,zb) … … 405 398 ! z1(ig) = pcapcal(ig)*ptsrf(ig) + cpp*zb(ig,1)*zc(ig,1) 406 399 ! & + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep 407 ! z2(ig) = pcapcal(ig) + cpp*zb(ig,1)*(1.-zd(ig,1)) 400 ! z2(ig) = pcapcal(ig) + cpp*zb(ig,1)*(1.-zd(ig,1)) 408 401 ! & +zdplanck(ig) 409 402 ! ztsrf2(ig) = z1(ig) / z2(ig) … … 432 425 ! Calculate vertical flux from the bottom to the first layer (dust) 433 426 ! ----------------------------------------------------------------- 434 do ig=1,ngrid 427 do ig=1,ngrid 435 428 rho(ig) = zb0(ig,1) /ptimestep 436 429 end do … … 440 433 ! Implicit inversion of q 441 434 ! ----------------------- 442 do iq=1,nq 443 444 if (iq.ne.igcm_h2o_vap) then435 do iq=1,nq 436 437 ! if (iq.ne.igcm_h2o_vap) then 445 438 446 439 DO ig=1,ngrid … … 448 441 zcq(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) 449 442 zdq(ig,nlay)=zb(ig,nlay)*z1(ig) 450 ENDDO 451 443 ENDDO 444 452 445 DO ilay=nlay-1,2,-1 453 446 DO ig=1,ngrid … … 460 453 ENDDO 461 454 462 ! if ((water).and.(iq.eq.iice)) then463 ! ! special case for water ice tracer: do not include464 ! ! h2o ice tracer from surface (which is set when handling465 ! ! h2o vapour case (see further down).466 ! ! zb(ig,1)=0 if iq ne ivap467 ! DO ig=1,ngrid468 ! z1(ig)=1./(za(ig,1)+469 ! & zb(ig,2)*(1.-zdq(ig,2)))470 ! zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+471 ! & zb(ig,2)*zcq(ig,2))*z1(ig)472 ! ENDDO473 ! else ! general case474 455 DO ig=1,ngrid 475 456 z1(ig)=1./(za(ig,1)+ … … 497 478 end do 498 479 499 endif ! if (iq.ne.igcm_h2o_vap)500 501 480 ! Removed (AF24): Calculate temperature tendency including latent heat term 502 481 ! and assuming an infinite source of water on the ground … … 517 496 pdvdif(ig,ilev)=(zv(ig,ilev)- 518 497 & (pv(ig,ilev)))/ptimestep 519 hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep 498 hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep 520 499 521 500 pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep … … 525 504 DO ig=1,ngrid ! computing sensible heat flux (atm => surface) 526 505 sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zh(ig,1)-ztsrf2(ig)) 527 ENDDO 506 ENDDO 528 507 529 508 if (tracer) then … … 561 540 ! endif 562 541 563 endif 542 endif 564 543 565 544 ! if(water)then
Note: See TracChangeset
for help on using the changeset viewer.