MODULE inifis_mod IMPLICIT NONE CONTAINS SUBROUTINE inifis(ngrid,nlayer,nq, & day_ini,pdaysec,nday,ptimestep, & plat,plon,parea, & prad,pg,pr,pcpp) use init_print_control_mod, only: init_print_control use radinc_h, only: ini_radinc_h, naerkind use datafile_mod, only: datadir use comdiurn_h, only: sinlat, coslat, sinlon, coslon use comgeomfi_h, only: totarea, totarea_planet use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, & init_time, daysec, dtphys use comcstfi_mod, only: rad, cpp, g, r, rcp, & mugaz, pi, avocado use planete_mod, only: nres use planetwide_mod, only: planetwide_sumval use callkeys_mod use ioipsl_getin_p_mod, only : getin_p use mod_phys_lmdz_para, only : is_parallel, is_master, bcast !======================================================================= ! ! purpose: ! ------- ! ! Initialisation for the physical parametrisations of the LMD ! Generic Model. ! ! author: Frederic Hourdin 15 / 10 /93 ! ------- ! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def) ! Ehouarn Millour (oct. 2008) tracers are now identified ! by their names and may not be contiguously ! stored in the q(:,:,:,:) array ! E.M. (june 2009) use getin routine to load parameters ! ! ! arguments: ! ---------- ! ! input: ! ------ ! ! ngrid Size of the horizontal grid. ! All internal loops are performed on that grid. ! nlayer Number of vertical layers. ! pdayref Day of reference for the simulation ! pday Number of days counted from the North. Spring ! equinoxe. ! !======================================================================= ! !----------------------------------------------------------------------- ! declarations: ! ------------- use datafile_mod, only: datadir use ioipsl_getin_p_mod, only: getin_p IMPLICIT NONE REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep INTEGER,INTENT(IN) :: nday INTEGER,INTENT(IN) :: ngrid,nlayer,nq REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid) integer,intent(in) :: day_ini INTEGER ig CHARACTER(len=20) :: rname="inifis" ! routine name, for messages EXTERNAL iniorbit,orbite EXTERNAL SSUM REAL SSUM ! Initialize flags lunout, prt_level, debug (in print_control_mod) CALL init_print_control ! initialize constants in comcstfi_mod rad=prad cpp=pcpp g=pg r=pr rcp=r/cpp mugaz=8.314*1000./pr ! dummy init pi=2.*asin(1.) avocado = 6.02214179e23 ! added by RW ! Initialize some "temporal and calendar" related variables #ifndef MESOSCALE CALL init_time(day_ini,pdaysec,nday,ptimestep) #endif ! read in some parameters from "run.def" for physics, ! or shared between dynamics and physics. call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics, ! in dynamical steps call getin_p("day_step",day_step) ! number of dynamical steps per day call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step ! do we read a startphy.nc file? (default: .true.) call getin_p("startphy_file",startphy_file) ! -------------------------------------------------------------- ! Reading the "callphys.def" file controlling some key options ! -------------------------------------------------------------- IF (is_master) THEN ! check if 'callphys.def' file is around inquire(file="callphys.def",exist=iscallphys) write(*,*) trim(rname)//": iscallphys=",iscallphys ENDIF call bcast(iscallphys) IF(iscallphys) THEN if (is_master) then write(*,*) write(*,*)'--------------------------------------------' write(*,*)trim(rname)//': Parameters for the physics (callphys.def)' write(*,*)'--------------------------------------------' endif if (is_master) write(*,*) trim(rname)//& ": Directory where external input files are:" ! default 'datadir' is set in "datadir_mod" call getin_p("datadir",datadir) ! default path if (is_master) write(*,*) trim(rname)//": datadir = ",trim(datadir) if (is_master) write(*,*) trim(rname)//& ": Run with or without tracer transport ?" tracer=.false. ! default value call getin_p("tracer",tracer) if (is_master) write(*,*) trim(rname)//": tracer = ",tracer if (is_master) write(*,*) trim(rname)//& ": Run with or without atm mass update "//& " due to tracer evaporation/condensation?" mass_redistrib=.false. ! default value call getin_p("mass_redistrib",mass_redistrib) if (is_master) write(*,*) trim(rname)//& ": mass_redistrib = ",mass_redistrib if (is_master) then write(*,*) trim(rname)//": Diurnal cycle ?" write(*,*) trim(rname)//& ": (if diurnal=false, diurnal averaged solar heating)" endif diurnal=.true. ! default value call getin_p("diurnal",diurnal) if (is_master) write(*,*) trim(rname)//": diurnal = ",diurnal if (is_master) then write(*,*) trim(rname)//": Seasonal cycle ?" write(*,*) trim(rname)//& ": (if season=false, Ls stays constant, to value ", & "set in 'start'" endif season=.true. ! default value call getin_p("season",season) if (is_master) write(*,*) trim(rname)//": season = ",season if (is_master) write(*,*) trim(rname)//& ": No seasonal cycle: initial day to lock the run during restart" noseason_day=0.0 ! default value call getin_p("noseason_day",noseason_day) if (is_master) write(*,*) trim(rname)//": noseason_day=", noseason_day if (is_master) write(*,*) trim(rname)//": Tidally resonant rotation ?" tlocked=.false. ! default value call getin_p("tlocked",tlocked) if (is_master) write(*,*) trim(rname)//": tlocked = ",tlocked if (is_master) write(*,*) trim(rname)//": Saturn ring shadowing ?" rings_shadow = .false. call getin_p("rings_shadow", rings_shadow) if (is_master) write(*,*) trim(rname)//": rings_shadow = ", rings_shadow if (is_master) write(*,*) trim(rname)//& ": Compute latitude-dependent gravity field?" oblate = .false. call getin_p("oblate", oblate) if (is_master) write(*,*) trim(rname)//": oblate = ", oblate if (is_master) write(*,*) trim(rname)//": Flattening of the planet (a-b)/a " flatten = 0.0 call getin_p("flatten", flatten) if (is_master) write(*,*) trim(rname)//": flatten = ", flatten if (is_master) write(*,*) trim(rname)//": Needed if oblate=.true.: J2" J2 = 0.0 call getin_p("J2", J2) if (is_master) write(*,*) trim(rname)//": J2 = ", J2 if (is_master) write(*,*) trim(rname)//& ": Needed if oblate=.true.: Planet mass (*1e24 kg)" MassPlanet = 0.0 call getin_p("MassPlanet", MassPlanet) if (is_master) write(*,*) trim(rname)//": MassPlanet = ", MassPlanet if (is_master) write(*,*) trim(rname)//& ": Needed if oblate=.true.: Planet mean radius (m)" Rmean = 0.0 call getin_p("Rmean", Rmean) if (is_master) write(*,*) trim(rname)//": Rmean = ", Rmean ! Test of incompatibility: ! if tlocked, then diurnal should be false if (tlocked.and.diurnal) then call abort_physic(rname,"If diurnal=true, we should turn off tlocked",1) endif if (is_master) write(*,*) trim(rname)//": Tidal resonance ratio ?" nres=0 ! default value call getin_p("nres",nres) if (is_master) write(*,*) trim(rname)//": nres = ",nres if (is_master) write(*,*) trim(rname)//& ": Write some extra output to the screen ?" lwrite=.false. ! default value call getin_p("lwrite",lwrite) if (is_master) write(*,*) trim(rname)//": lwrite = ",lwrite if (is_master) write(*,*) trim(rname)//& ": Save statistics in file stats.nc ?" callstats=.true. ! default value call getin_p("callstats",callstats) if (is_master) write(*,*) trim(rname)//": callstats = ",callstats if (is_master) write(*,*) trim(rname)//& ": Test energy conservation of model physics ?" enertest=.false. ! default value call getin_p("enertest",enertest) if (is_master) write(*,*) trim(rname)//": enertest = ",enertest if (is_master) write(*,*) trim(rname)//": call radiative transfer ?" callrad=.true. ! default value call getin_p("callrad",callrad) if (is_master) write(*,*) trim(rname)//": callrad = ",callrad if (is_master) write(*,*) trim(rname)//& ": call correlated-k radiative transfer ?" corrk=.true. ! default value call getin_p("corrk",corrk) if (is_master) write(*,*) trim(rname)//": corrk = ",corrk if (is_master) write(*,*) trim(rname)//& ": prohibit calculations outside corrk T grid?" strictboundcorrk=.true. ! default value call getin_p("strictboundcorrk",strictboundcorrk) if (is_master) write(*,*) trim(rname)//& ": strictboundcorrk = ",strictboundcorrk if (is_master) write(*,*) trim(rname)//& ": prohibit calculations outside CIA T grid?" strictboundcia=.true. ! default value call getin_p("strictboundcia",strictboundcia) if (is_master) write(*,*) trim(rname)//& ": strictboundcia = ",strictboundcia if (is_master) write(*,*) trim(rname)//& ": Minimum atmospheric temperature for Planck function integration ?" tplanckmin=30.0 ! default value call getin_p("tplanckmin",tplanckmin) if (is_master) write(*,*) trim(rname)//": tplanckmin = ",tplanckmin if (is_master) write(*,*) trim(rname)//& ": Maximum atmospheric temperature for Planck function integration ?" tplanckmax=1500.0 ! default value call getin_p("tplanckmax",tplanckmax) if (is_master) write(*,*) trim(rname)//": tplanckmax = ",tplanckmax if (is_master) write(*,*) trim(rname)//& ": Temperature step for Planck function integration ?" dtplanck=0.1 ! default value call getin_p("dtplanck",dtplanck) if (is_master) write(*,*) trim(rname)//": dtplanck = ",dtplanck if (is_master) write(*,*) trim(rname)//& ": call gaseous absorption in the visible bands?"//& " (matters only if callrad=T)" callgasvis=.false. ! default value call getin_p("callgasvis",callgasvis) if (is_master) write(*,*) trim(rname)//": callgasvis = ",callgasvis if (is_master) write(*,*) trim(rname)//& ": call continuum opacities in radiative transfer ?"//& " (matters only if callrad=T)" continuum=.true. ! default value call getin_p("continuum",continuum) if (is_master) write(*,*) trim(rname)//": continuum = ",continuum if (is_master) write(*,*) trim(rname)//": version for H2H2 CIA file ?" versH2H2cia=2011 ! default value (should be 2018 but retrocompatibility first) call getin_p("versH2H2cia",versH2H2cia) if (is_master) write(*,*) trim(rname)//": versH2H2cia = ",versH2H2cia ! Sanity check if (versH2H2cia.ne.2011 .and. versH2H2cia.ne.2018) then call abort_physic(rname,"Error: Choose a correct value (2011 or 2018) for versH2H2cia !",1) endif if (is_master) write(*,*) trim(rname)//& ": CIA - normal or equilibrium ortho-para mixture for H2?" H2orthopara_mixture = 'normal' call getin_p("H2orthopara_mixture",H2orthopara_mixture) if (is_master) write(*,*)trim(rname)//& ": H2orthopara_mixture = ",trim(H2orthopara_mixture) if (is_master) write(*,*) trim(rname)//& ": call turbulent vertical diffusion ?" calldifv=.true. ! default value call getin_p("calldifv",calldifv) if (is_master) write(*,*) trim(rname)//": calldifv = ",calldifv if (is_master) write(*,*) trim(rname)//": use turbdiff instead of vdifc ?" UseTurbDiff=.true. ! default value call getin_p("UseTurbDiff",UseTurbDiff) if (is_master) write(*,*) trim(rname)//": UseTurbDiff = ",UseTurbDiff if (is_master) write(*,*) trim(rname)//": call convective adjustment ?" calladj=.true. ! default value call getin_p("calladj",calladj) if (is_master) write(*,*) trim(rname)//": calladj = ",calladj if (is_master) write(*,*)trim(rname)//& ": call Lott's non-oro GWs parameterisation scheme ?" calllott_nonoro=.false. ! default value call getin_p("calllott_nonoro",calllott_nonoro) if (is_master) write(*,*)trim(rname)//& ": calllott_nonoro = ",calllott_nonoro if (is_master) write(*,*) trim(rname)//": call thermal plume model ?" calltherm=.false. ! default value call getin_p("calltherm",calltherm) if (is_master) write(*,*) trim(rname)//": calltherm = ",calltherm if (is_master) write(*,*) trim(rname)//": call CO2 condensation ?" co2cond=.false. ! default value call getin_p("co2cond",co2cond) if (is_master) write(*,*) trim(rname)//": co2cond = ",co2cond ! Test of incompatibility if (co2cond.and.(.not.tracer)) then call abort_physic(rname,"Error: We need a CO2 ice tracer to condense CO2!",1) endif if (is_master) write(*,*) trim(rname)//": CO2 supersaturation level ?" co2supsat=1.0 ! default value call getin_p("co2supsat",co2supsat) if (is_master) write(*,*) trim(rname)//": co2supsat = ",co2supsat if (is_master) write(*,*) trim(rname)//& ": Radiative timescale for Newtonian cooling (in Earth days)?" tau_relax=30. ! default value call getin_p("tau_relax",tau_relax) if (is_master) write(*,*) trim(rname)//": tau_relax = ",tau_relax tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds if (is_master) write(*,*)trim(rname)//& ": call thermal conduction in the soil ?" callsoil=.true. ! default value call getin_p("callsoil",callsoil) if (is_master) write(*,*) trim(rname)//": callsoil = ",callsoil if (is_master) write(*,*)trim(rname)//& ": Rad transfer is computed every iradia", & " physical timestep" iradia=1 ! default value call getin_p("iradia",iradia) if (is_master) write(*,*)trim(rname)//": iradia = ",iradia if (is_master) write(*,*)trim(rname)//": Rayleigh scattering ?" rayleigh=.false. call getin_p("rayleigh",rayleigh) if (is_master) write(*,*)trim(rname)//": rayleigh = ",rayleigh if (is_master) write(*,*) trim(rname)//& ": Use blackbody for stellar spectrum ?" stelbbody=.false. ! default value call getin_p("stelbbody",stelbbody) if (is_master) write(*,*) trim(rname)//": stelbbody = ",stelbbody if (is_master) write(*,*) trim(rname)//": Stellar blackbody temperature ?" stelTbb=5800.0 ! default value call getin_p("stelTbb",stelTbb) if (is_master) write(*,*) trim(rname)//": stelTbb = ",stelTbb if (is_master) write(*,*) trim(rname)//": Output mean OLR in 1D?" meanOLR=.false. call getin_p("meanOLR",meanOLR) if (is_master) write(*,*) trim(rname)//": meanOLR = ",meanOLR if (is_master) write(*,*)trim(rname)//": Output spectral OLR in 3D?" specOLR=.false. call getin_p("specOLR",specOLR) if (is_master) write(*,*)trim(rname)//": specOLR = ",specOLR if (is_master) write(*,*)trim(rname)//& ": Output diagnostic optical thickness attenuation exp(-tau)?" diagdtau=.false. call getin_p("diagdtau",diagdtau) if (is_master) write(*,*)trim(rname)//": diagdtau = ",diagdtau if (is_master) write(*,*)trim(rname)//": Operate in kastprof mode?" kastprof=.false. call getin_p("kastprof",kastprof) if (is_master) write(*,*)trim(rname)//": kastprof = ",kastprof if (is_master) write(*,*)trim(rname)//& ": Uniform absorption in radiative transfer?" graybody=.false. call getin_p("graybody",graybody) if (is_master) write(*,*)trim(rname)//": graybody = ",graybody ! Soil model if (is_master) write(*,*)trim(rname)//& ": Number of sub-surface layers for soil scheme?" ! default value of nsoilmx set in comsoil_h call getin_p("nsoilmx",nsoilmx) if (is_master) write(*,*)trim(rname)//": nsoilmx=",nsoilmx if (is_master) write(*,*)trim(rname)//& ": Thickness of topmost soil layer (m)?" ! default value of lay1_soil set in comsoil_h call getin_p("lay1_soil",lay1_soil) if (is_master) write(*,*)trim(rname)//": lay1_soil = ",lay1_soil if (is_master) write(*,*)trim(rname)//& ": Coefficient for soil layer thickness distribution?" ! default value of alpha_soil set in comsoil_h call getin_p("alpha_soil",alpha_soil) if (is_master) write(*,*)trim(rname)//": alpha_soil = ",alpha_soil ! Slab Ocean if (is_master) write(*,*) trim(rname)//": Use slab-ocean ?" ok_slab_ocean=.false. ! default value call getin_p("ok_slab_ocean",ok_slab_ocean) if (is_master) write(*,*) trim(rname)//": ok_slab_ocean = ",ok_slab_ocean ! Sanity check: for now slab oncean only works in serial mode if (ok_slab_ocean.and.is_parallel) then call abort_physic(rname," Error: slab ocean should only be used in serial mode!",1) endif if (is_master) write(*,*) trim(rname)//": Use slab-sea-ice ?" ok_slab_sic=.true. ! default value call getin_p("ok_slab_sic",ok_slab_sic) if (is_master) write(*,*) trim(rname)//": ok_slab_sic = ",ok_slab_sic if (is_master) write(*,*) trim(rname)//& ": Use heat transport for the ocean ?" ok_slab_heat_transp=.true. ! default value call getin_p("ok_slab_heat_transp",ok_slab_heat_transp) if (is_master) write(*,*) trim(rname)//& ": ok_slab_heat_transp = ",ok_slab_heat_transp ! Photochemistry and chemistry in the thermosphere if (is_master) write(*,*) trim(rname)//": Use photochemistry ?" photochem=.false. ! default value call getin_p("photochem",photochem) if (is_master) write(*,*) trim(rname)//": photochem = ",photochem if (is_master) write(*,*) trim(rname)//": Use photolysis heat table ?" photoheat=.false. ! default value call getin_p("photoheat",photoheat) if (is_master) write(*,*) trim(rname)//": photoheat = ",photoheat if (is_master) write(*,*) trim(rname)//& ": Use photolysis online calculation ?" jonline=.false. ! default value call getin_p("jonline",jonline) if (is_master) write(*,*) trim(rname)//": jonline = ",jonline if (is_master) write(*,*) trim(rname)//": Use deposition ?" depos=.false. ! default value call getin_p("depos",depos) if (is_master) write(*,*) trim(rname)//": depos = ",depos if (is_master) write(*,*)trim(rname)//": Production of haze ?" haze=.false. ! default value call getin_p("haze",haze) if (is_master) write(*,*)trim(rname)//": haze = ",haze ! Global1D mean and solar zenith angle if (ngrid.eq.1) then PRINT*, 'Simulate global averaged conditions ?' global1d = .false. ! default value call getin_p("global1d",global1d) write(*,*) "global1d = ",global1d ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle. if (global1d.and.diurnal) then call abort_physic(rname,'if global1d is true, diurnal must be set to false',1) endif if (global1d) then PRINT *,'Solar Zenith angle (deg.) ?' PRINT *,'(assumed for averaged solar flux S/4)' szangle=60.0 ! default value call getin_p("szangle",szangle) write(*,*) "szangle = ",szangle endif endif ! of if (ngrid.eq.1) ! Test of incompatibility: ! if kastprof used, we must be in 1D if (kastprof.and.(ngrid.gt.1)) then call abort_physic(rname,'kastprof can only be used in 1D!',1) endif if (is_master) write(*,*)trim(rname)//& ": Stratospheric temperature for kastprof mode?" Tstrat=167.0 call getin_p("Tstrat",Tstrat) if (is_master) write(*,*)trim(rname)//": Tstrat = ",Tstrat if (is_master) write(*,*)trim(rname)//": Remove lower boundary?" nosurf=.false. call getin_p("nosurf",nosurf) if (is_master) write(*,*)trim(rname)//": nosurf = ",nosurf ! Tests of incompatibility: if (nosurf.and.callsoil) then if (is_master) then write(*,*)trim(rname)//'nosurf not compatible with soil scheme!' write(*,*)trim(rname)//'... got to make a choice!' endif call abort_physic(rname,"nosurf not compatible with soil scheme!",1) endif if (is_master) write(*,*)trim(rname)//": Add an internal heat flux?", & "... matters only if callsoil=F" intheat=0. call getin_p("intheat",intheat) if (is_master) write(*,*)trim(rname)//": intheat = ",intheat if (is_master) write(*,*)trim(rname)//& ": Use Newtonian cooling for radiative transfer?" newtonian=.false. call getin_p("newtonian",newtonian) if (is_master) write(*,*)trim(rname)//": newtonian = ",newtonian ! Tests of incompatibility: if (newtonian.and.corrk) then call abort_physic(rname,'newtonian not compatible with correlated-k!',1) endif if (newtonian.and.calladj) then call abort_physic(rname,'newtonian not compatible with adjustment!',1) endif if (newtonian.and.calldifv) then call abort_physic(rname,'newtonian not compatible with a boundary layer!',1) endif if (is_master) write(*,*)trim(rname)//": Test physics timescale in 1D?" testradtimes=.false. call getin_p("testradtimes",testradtimes) if (is_master) write(*,*)trim(rname)//": testradtimes = ",testradtimes ! Test of incompatibility: ! if testradtimes used, we must be in 1D if (testradtimes.and.(ngrid.gt.1)) then call abort_physic(rname,'testradtimes can only be used in 1D!',1) endif if (is_master) write(*,*)trim(rname)//": Default planetary temperature?" tplanet=215.0 call getin_p("tplanet",tplanet) if (is_master) write(*,*)trim(rname)//": tplanet = ",tplanet if (is_master) write(*,*)trim(rname)//": Which star?" startype=1 ! default value = Sol call getin_p("startype",startype) if (is_master) write(*,*)trim(rname)//": startype = ",startype if (is_master) write(*,*)trim(rname)//": Value of stellar flux at 1 AU?" Fat1AU=1356.0 ! default value = Sol today call getin_p("Fat1AU",Fat1AU) if (is_master) write(*,*)trim(rname)//": Fat1AU = ",Fat1AU ! TRACERS: if (is_master) write(*,*)trim(rname)//": Varying H2O cloud fraction?" CLFvarying=.false. ! default value call getin_p("CLFvarying",CLFvarying) if (is_master) write(*,*)trim(rname)//": CLFvarying = ",CLFvarying if (is_master) write(*,*)trim(rname)//& ": Value of fixed H2O cloud fraction?" CLFfixval=1.0 ! default value call getin_p("CLFfixval",CLFfixval) if (is_master) write(*,*)trim(rname)//": CLFfixval = ",CLFfixval if (is_master) write(*,*)trim(rname)//": fixed radii for Cloud particles?" radfixed=.false. ! default value call getin_p("radfixed",radfixed) if (is_master) write(*,*)trim(rname)//": radfixed = ",radfixed if(kastprof)then radfixed=.true. endif if (is_master) write(*,*)trim(rname)//& ": Number mixing ratio of CO2 ice particles:" Nmix_co2=1.e6 ! default value call getin_p("Nmix_co2",Nmix_co2) if (is_master) write(*,*)trim(rname)//": Nmix_co2 = ",Nmix_co2 ! write(*,*)"Number of radiatively active aerosols:" ! naerkind=0. ! default value ! call getin_p("naerkind",naerkind) ! write(*,*)" naerkind = ",naerkind if (is_master) write(*,*)trim(rname)//": Opacity of dust (if used):" dusttau=0. ! default value call getin_p("dusttau",dusttau) if (is_master) write(*,*)trim(rname)//": dusttau = ",dusttau if (is_master) write(*,*)trim(rname)//": Radiatively active CO2 aerosols?" aeroco2=.false. ! default value call getin_p("aeroco2",aeroco2) if (is_master) write(*,*)trim(rname)//": aeroco2 = ",aeroco2 if (is_master) write(*,*)trim(rname)//": Fixed CO2 aerosol distribution?" aerofixco2=.false. ! default value call getin_p("aerofixco2",aerofixco2) if (is_master) write(*,*)trim(rname)//": aerofixco2 = ",aerofixco2 if (is_master) write(*,*)trim(rname)//": Radiatively active water ice?" aeroh2o=.false. ! default value call getin_p("aeroh2o",aeroh2o) if (is_master) write(*,*)trim(rname)//": aeroh2o = ",aeroh2o if (is_master) write(*,*)trim(rname)//": Fixed H2O aerosol distribution?" aerofixh2o=.false. ! default value call getin_p("aerofixh2o",aerofixh2o) if (is_master) write(*,*)trim(rname)//": aerofixh2o = ",aerofixh2o if (is_master) write(*,*)trim(rname)//& ": Radiatively active sulfuric acid aerosols?" aeroh2so4=.false. ! default value call getin_p("aeroh2so4",aeroh2so4) if (is_master) write(*,*)trim(rname)//": aeroh2so4 = ",aeroh2so4 if (is_master) write(*,*)trim(rname)//& ": Radiatively active auroral aerosols?" aeroaurora=.false. ! default value call getin_p("aeroaurora",aeroaurora) if (is_master) write(*,*)trim(rname)//": aeroaurora = ",aeroaurora !================================= ! TWOLAY scheme and NH3 cloudare left for retrocompatibility only, ! You should now use N-LAYER scheme (see below). if (is_master) write(*,*)trim(rname)//& ": Radiatively active two-layer aerosols?" aeroback2lay=.false. ! default value call getin_p("aeroback2lay",aeroback2lay) if (is_master) write(*,*)trim(rname)//": aeroback2lay = ",aeroback2lay if (aeroback2lay.and.is_master) then print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...' print*,'You should use the generic n-layer scheme (see aeronlay).' endif if (is_master) write(*,*)trim(rname)//": Radiatively active ammonia cloud?" aeronh3=.false. ! default value call getin_p("aeronh3",aeronh3) if (is_master) write(*,*)trim(rname)//": aeronh3 = ",aeronh3 if (aeronh3.and.is_master) then print*,'Warning : You are using specific NH3 cloud scheme ...' print*,'You should use the generic n-layer scheme (see aeronlay).' endif if (is_master) write(*,*)trim(rname)//& ": TWOLAY AEROSOL: total optical depth "//& "in the tropospheric layer (visible)" obs_tau_col_tropo=8.D0 call getin_p("obs_tau_col_tropo",obs_tau_col_tropo) if (is_master) write(*,*)trim(rname)//& ": obs_tau_col_tropo = ",obs_tau_col_tropo if (is_master) write(*,*)trim(rname)//& ": TWOLAY AEROSOL: total optical depth "//& "in the stratospheric layer (visible)" obs_tau_col_strato=0.08D0 call getin_p("obs_tau_col_strato",obs_tau_col_strato) if (is_master) write(*,*)trim(rname)//& ": obs_tau_col_strato = ",obs_tau_col_strato if (is_master) write(*,*)trim(rname)//& ": TWOLAY AEROSOL: optprop_back2lay_vis?" optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat' call getin_p("optprop_back2lay_vis",optprop_back2lay_vis) if (is_master) write(*,*)trim(rname)//& ": optprop_back2lay_vis = ",trim(optprop_back2lay_vis) if (is_master) write(*,*)trim(rname)//& ": TWOLAY AEROSOL: optprop_back2lay_ir?" optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat' call getin_p("optprop_back2lay_ir",optprop_back2lay_ir) if (is_master) write(*,*)trim(rname)//& ": optprop_back2lay_ir = ",trim(optprop_back2lay_ir) if (is_master) write(*,*)trim(rname)//& ": TWOLAY AEROSOL: pres_bottom_tropo? in pa" pres_bottom_tropo=66000.0 call getin_p("pres_bottom_tropo",pres_bottom_tropo) if (is_master) write(*,*)trim(rname)//& ": pres_bottom_tropo = ",pres_bottom_tropo if (is_master) write(*,*)trim(rname)//& ": TWOLAY AEROSOL: pres_top_tropo? in pa" pres_top_tropo=18000.0 call getin_p("pres_top_tropo",pres_top_tropo) if (is_master) write(*,*)trim(rname)//& ": pres_top_tropo = ",pres_top_tropo if (is_master) write(*,*)trim(rname)//& ": TWOLAY AEROSOL: pres_bottom_strato? in pa" pres_bottom_strato=2000.0 call getin_p("pres_bottom_strato",pres_bottom_strato) if (is_master) write(*,*)trim(rname)//& ": pres_bottom_strato = ",pres_bottom_strato ! Sanity check if (pres_bottom_strato .gt. pres_top_tropo) then if(is_master) then print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def' print*,'you have pres_top_tropo > pres_bottom_strato !' endif call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1) endif if (is_master) write(*,*)trim(rname)//& ": TWOLAY AEROSOL: pres_top_strato? in pa" pres_top_strato=100.0 call getin_p("pres_top_strato",pres_top_strato) if (is_master) write(*,*)trim(rname)//& ": pres_top_strato = ",pres_top_strato if (is_master) write(*,*)trim(rname)//& ": TWOLAY AEROSOL: particle size in the ", & "tropospheric layer, in meters" size_tropo=2.e-6 call getin_p("size_tropo",size_tropo) if (is_master) write(*,*)trim(rname)//": size_tropo = ",size_tropo if (is_master) write(*,*)trim(rname)//& ": TWOLAY AEROSOL: particle size in the ", & "stratospheric layer, in meters" size_strato=1.e-7 call getin_p("size_strato",size_strato) if (is_master) write(*,*)trim(rname)//": size_strato = ",size_strato if (is_master) write(*,*)trim(rname)//& ": NH3 (thin) cloud: total optical depth" tau_nh3_cloud=7.D0 call getin_p("tau_nh3_cloud",tau_nh3_cloud) if (is_master) write(*,*)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud if (is_master) write(*,*)trim(rname)//": NH3 (thin) cloud pressure level" pres_nh3_cloud=7.D0 call getin_p("pres_nh3_cloud",pres_nh3_cloud) if (is_master) write(*,*)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud if (is_master) write(*,*)trim(rname)//": NH3 (thin) cloud: particle sizes" size_nh3_cloud=3.e-6 call getin_p("size_nh3_cloud",size_nh3_cloud) if (is_master) write(*,*)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud !================================= ! Generic N-LAYER aerosol scheme if (is_master) write(*,*)trim(rname)//& ": Radiatively active generic n-layer aerosols?" aeronlay=.false. ! default value call getin_p("aeronlay",aeronlay) if (is_master) write(*,*)trim(rname)//": aeronlay = ",aeronlay if (is_master) write(*,*)trim(rname)//& ": Number of generic aerosols layers?" nlayaero=1 ! default value call getin_p("nlayaero",nlayaero) ! Avoid to allocate arrays of size 0 if (aeronlay .and. nlayaero.lt.1) then if (is_master) then print*, " You are trying to set no generic aerosols..." print*, " Set aeronlay=.false. instead ! I abort." endif call abort_physic(rname,"no generic aerosols...",1) endif if (.not. aeronlay) nlayaero=1 if (is_master) write(*,*)trim(rname)//": nlayaero = ",nlayaero ! This is necessary, we just set the number of aerosol layers IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero)) IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero)) IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero)) IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero)) IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero)) IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero)) IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero)) IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero)) IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero)) IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero)) if (is_master) write(*,*)trim(rname)//& ": Generic n-layer aerosols: Optical depth at reference wavelenght" aeronlay_tauref=1.0E-1 call getin_p("aeronlay_tauref",aeronlay_tauref) if (is_master) write(*,*)trim(rname)//& ": aeronlay_tauref = ",aeronlay_tauref if (is_master) write(*,*)trim(rname)//& ": Generic n-layer aerosols: Reference wavelenght for optical depths (m)" aeronlay_lamref=0.6E-6 call getin_p("aeronlay_lamref",aeronlay_lamref) if (is_master) write(*,*)trim(rname)//& ": aeronlay_lamref = ",aeronlay_lamref if (is_master) then write(*,*)trim(rname)//& ": Generic n-layer aerosols: Vertical profile choice : " write(*,*)trim(rname)//& " (1) Tau btwn ptop and pbot follows atm. scale height" write(*,*)trim(rname)//& " or (2) Tau above pbot follows its own scale height" endif aeronlay_choice=1 call getin_p("aeronlay_choice",aeronlay_choice) if (is_master) write(*,*)trim(rname)//& ": aeronlay_choice = ",aeronlay_choice if (is_master) write(*,*)trim(rname)//& ": Generic n-layer aerosols: bottom pressures (Pa)" aeronlay_pbot=2000.0 call getin_p("aeronlay_pbot",aeronlay_pbot) if (is_master) write(*,*)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot if (is_master) write(*,*)trim(rname)//& ": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) " aeronlay_ptop=300000.0 call getin_p("aeronlay_ptop",aeronlay_ptop) if (is_master) write(*,*)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop if (is_master) write(*,*)trim(rname)//& ": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height" aeronlay_sclhght=0.2 call getin_p("aeronlay_sclhght",aeronlay_sclhght) if (is_master) write(*,*)trim(rname)//& ": aeronlay_sclhght = ",aeronlay_sclhght if (is_master) write(*,*)trim(rname)//& ": Generic n-layer aerosols: particles effective radii (m)" aeronlay_size=1.e-6 call getin_p("aeronlay_size",aeronlay_size) if (is_master) write(*,*)trim(rname)//": aeronlay_size = ",aeronlay_size if (is_master) write(*,*)trim(rname)//& ": Generic n-layer aerosols: particles radii effective variance" aeronlay_nueff=0.1 call getin_p("aeronlay_nueff",aeronlay_nueff) if (is_master) write(*,*)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff if (is_master) write(*,*)trim(rname)//& ": Generic n-layer aerosols: VIS optical properties file" optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat' call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis) if (is_master) write(*,*)trim(rname)//& ": optprop_aeronlay_vis = ",optprop_aeronlay_vis if (is_master) write(*,*)trim(rname)//& ": Generic n-layer aerosols: IR optical properties file" optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat' call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir) if (is_master) write(*,*)trim(rname)//& ": optprop_aeronlay_ir = ",optprop_aeronlay_ir !================================= if (is_master) write(*,*)trim(rname)//& ": Cloud pressure level (with kastprof only):" cloudlvl=0. ! default value call getin_p("cloudlvl",cloudlvl) if (is_master) write(*,*)trim(rname)//": cloudlvl = ",cloudlvl if (is_master) write(*,*)trim(rname)//& ": Is the variable gas species radiatively active?" Tstrat=167.0 varactive=.false. call getin_p("varactive",varactive) if (is_master) write(*,*)trim(rname)//": varactive = ",varactive if (is_master) write(*,*)trim(rname)//& ": Is the variable gas species distribution set?" varfixed=.false. call getin_p("varfixed",varfixed) if (is_master) write(*,*)trim(rname)//": varfixed = ",varfixed if (is_master) write(*,*)trim(rname)//& ": What is the saturation % of the variable species?" satval=0.8 call getin_p("satval",satval) if (is_master) write(*,*)trim(rname)//": satval = ",satval ! Test of incompatibility: ! if varactive, then varfixed should be false if (varactive.and.varfixed) then call abort_physic(rname,'if varactive, varfixed must be OFF!',1) endif if (is_master) write(*,*)trim(rname)//": Gravitationnal sedimentation ?" sedimentation=.false. ! default value call getin_p("sedimentation",sedimentation) if (is_master) write(*,*)trim(rname)//": sedimentation = ",sedimentation if (is_master) write(*,*)trim(rname)//": Compute water cycle ?" water=.false. ! default value call getin_p("water",water) if (is_master) write(*,*)trim(rname)//": water = ",water ! Test of incompatibility: ! if water is true, there should be at least a tracer if (water.and.(.not.tracer)) then call abort_physic(rname,'if water is ON, tracer must be ON too!',1) endif if (is_master) write(*,*)trim(rname)//": Include water condensation ?" watercond=.false. ! default value call getin_p("watercond",watercond) if (is_master) write(*,*)trim(rname)//": watercond = ",watercond ! Test of incompatibility: ! if watercond is used, then water should be used too if (watercond.and.(.not.water)) then call abort_physic(rname,& 'if watercond is used, water should be used too',1) endif if (is_master) write(*,*)trim(rname)//": Include water precipitation ?" waterrain=.false. ! default value call getin_p("waterrain",waterrain) if (is_master) write(*,*)trim(rname)//": waterrain = ",waterrain if (is_master) write(*,*)trim(rname)//": Include surface hydrology ?" hydrology=.false. ! default value call getin_p("hydrology",hydrology) if (is_master) write(*,*)trim(rname)//": hydrology = ",hydrology if (is_master) write(*,*)trim(rname)//": Evolve surface water sources ?" sourceevol=.false. ! default value call getin_p("sourceevol",sourceevol) if (is_master) write(*,*)trim(rname)//": sourceevol = ",sourceevol if (is_master) write(*,*)trim(rname)//": Ice evolution timestep ?" icetstep=100.0 ! default value call getin_p("icetstep",icetstep) if (is_master) write(*,*)trim(rname)//": icetstep = ",icetstep if (is_master) write(*,*)trim(rname)//": Spectral Dependant albedo ?" albedo_spectral_mode=.false. ! default value call getin_p("albedo_spectral_mode",albedo_spectral_mode) if (is_master) write(*,*)trim(rname)//& ": albedo_spectral_mode = ",albedo_spectral_mode if (is_master) then write(*,*)trim(rname)//": Snow albedo ?" write(*,*)trim(rname)//": If albedo_spectral_mode=.true., then this " write(*,*)trim(rname)//": corresponds to the 0.5 microns snow albedo." endif albedosnow=0.5 ! default value call getin_p("albedosnow",albedosnow) if (is_master) write(*,*)trim(rname)//": albedosnow = ",albedosnow if (is_master) write(*,*)trim(rname)//": Ocean albedo ?" alb_ocean=0.07 ! default value call getin_p("alb_ocean",alb_ocean) if (is_master) write(*,*)trim(rname)//": alb_ocean = ",alb_ocean if (is_master) write(*,*)trim(rname)//": CO2 ice albedo ?" albedoco2ice=0.5 ! default value call getin_p("albedoco2ice",albedoco2ice) if (is_master) write(*,*)trim(rname)//": albedoco2ice = ",albedoco2ice if (is_master) write(*,*)trim(rname)//": Maximum ice thickness ?" maxicethick=2.0 ! default value call getin_p("maxicethick",maxicethick) if (is_master) write(*,*)trim(rname)//": maxicethick = ",maxicethick if (is_master) write(*,*)trim(rname)//": Freezing point of seawater ?" Tsaldiff=-1.8 ! default value call getin_p("Tsaldiff",Tsaldiff) if (is_master) write(*,*)trim(rname)//": Tsaldiff = ",Tsaldiff if (is_master) write(*,*)trim(rname)//": Minimum eddy mix coeff in 1D ?" kmixmin=1.0e-2 ! default value call getin_p("kmixmin",kmixmin) if (is_master) write(*,*)trim(rname)//": kmixmin = ",kmixmin if (is_master) write(*,*)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1' if (is_master) write(*,*)'Predefined Mg from dynamics is ',mugaz,'amu' force_cpp=.false. ! default value call getin_p("force_cpp",force_cpp) if (force_cpp) then if (is_master) write(*,*)trim(rname)//": force_cpp = ",force_cpp if (is_master) write(*,*)trim(rname)//": force_cpp is deprecated.",& "Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true." call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1) endif if (is_master) write(*,*)trim(rname)//& ": where do you want your cpp/mugaz value to come from?",& "=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?" cpp_mugaz_mode = 0 ! default value call getin_p("cpp_mugaz_mode",cpp_mugaz_mode) if (is_master) write(*,*)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then write(*,*)' !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics' write(*,*)' in 3D can result in very pathological behavior. You have been warned !!!!!' else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then ! for this specific 1D error we will remove run.def before aborting JL22 call system("rm -rf run.def") call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1) endif if (cpp_mugaz_mode == 1) then mugaz = -99999. if (is_master) write(*,*)trim(rname)//& ": MEAN MOLECULAR MASS in g mol-1 ?" call getin_p("mugaz",mugaz) IF (mugaz.eq.-99999.) THEN call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1) ENDIF cpp = -99999. if (is_master) write(*,*)trim(rname)//& ": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?" call getin_p("cpp",cpp) IF (cpp.eq.-99999.) THEN PRINT *, "cpp must be set if cpp_mugaz_mode = 1" STOP ENDIF if (is_master) write(*,*)'New Cp from callfis is ',cpp,'J kg^-1 K^-1' if (is_master) write(*,*)'New Mg from callfis is ',mugaz,'amu' endif ! of if (cpp_mugaz_mode == 1) call su_gases call calc_cpp_mugaz if (is_master) then PRINT*,'--------------------------------------------' PRINT* PRINT* endif ELSE call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1) ENDIF ! of IF(iscallphys) if (is_master) then PRINT* PRINT*,'inifis: daysec',daysec PRINT* PRINT*,'inifis: The radiative transfer is computed:' PRINT*,' each ',iradia,' physical time-step' PRINT*,' or each ',iradia*dtphys,' seconds' PRINT* endif !----------------------------------------------------------------------- ! Some more initialization: ! ------------------------ ! Initializations for comgeomfi_h #ifndef MESOSCALE totarea=SSUM(ngrid,parea,1) call planetwide_sumval(parea,totarea_planet) !! those are defined in comdiurn_h.F90 IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid)) IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid)) IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid)) IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid)) DO ig=1,ngrid sinlat(ig)=sin(plat(ig)) coslat(ig)=cos(plat(ig)) sinlon(ig)=sin(plon(ig)) coslon(ig)=cos(plon(ig)) ENDDO #endif ! initialize variables in radinc_h call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck) ! allocate "comsoil_h" arrays call ini_comsoil_h(ngrid) END SUBROUTINE inifis END MODULE inifis_mod