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 radcommon_h, only: ini_radcommon_h use radii_mod, only: radfixed, Nmix_n2 use datafile_mod, only: datadir,hazeprop_file,hazerad_file,hazemmr_file,hazedens_file 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 surfdat_h use wstats_mod, only: callstats 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(*,*) "Haze optical properties datafile" hazeprop_file="optprop_tholins_fractal050nm.dat" ! default file call getin_p("hazeprop_file",hazeprop_file) if (is_master) write(*,*) trim(rname)//" hazeprop_file = ",trim(hazeprop_file) if (is_master) write(*,*) "Haze radii datafile" hazerad_file="hazerad.txt" ! default file call getin_p("hazerad_file",hazerad_file) if (is_master) write(*,*) trim(rname)//" hazerad_file = ",trim(hazerad_file) if (is_master) write(*,*) "Haze mmr datafile" hazemmr_file="None" ! default file call getin_p("hazemmr_file",hazemmr_file) if (is_master) write(*,*) trim(rname)//" hazemmr_file = ",trim(hazemmr_file) if (is_master) write(*,*) "Haze density datafile" hazedens_file="None" ! default file call getin_p("hazedens_file",hazedens_file) if (is_master) write(*,*) trim(rname)//" hazedens_file = ",trim(hazedens_file) 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) then write(*,*) trim(rname)//": Fast mode (nogcm) ?" endif fast=.false. ! default value call getin_p("fast",fast) if (is_master) write(*,*) trim(rname)//": fast = ",fast if (is_master) write(*,*) trim(rname)//"Run with Triton orbit ?" triton=.false. ! default value call getin_p("triton",triton) if (is_master) write(*,*) trim(rname)//" triton = ",triton convergeps=.false. ! default value call getin_p("convergeps",convergeps) if (is_master) write(*,*) trim(rname)//" convergeps = ",convergeps conservn2=.false. ! default value call getin_p("conservn2",conservn2) if (is_master) write(*,*) trim(rname)//" conservn2 = ",conservn2 if (is_master) write(*,*) trim(rname)//"KBO runs (eris, makemake) ?" kbo=.false. ! default value call getin_p("kbo",kbo) if (is_master) write(*,*) trim(rname)//" kbo = ",kbo if (is_master) write(*,*) trim(rname)//"Specific paleo run ?" paleo=.false. ! default value call getin_p("paleo",paleo) if (is_master) write(*,*) trim(rname)//" paleo = ",paleo if (is_master) write(*,*) trim(rname)//"paleoclimate step (Earth years) " paleoyears=10000. ! default value call getin_p("paleoyears",paleoyears) if (is_master) write(*,*) trim(rname)//" paleoyears = ",paleoyears 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)//": 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)//& ": 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)//& ": 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 ! Chemistry in the thermosphere 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 if (is_master) write(*,*)trim(rname)// "call thermal conduction ?" callconduct=.false. ! default value call getin_p("callconduct",callconduct) if (is_master) write(*,*)trim(rname)// " callconduct = ",callconduct if (is_master) write(*,*)trim(rname)// "call phitop ?" phitop=0. ! default value call getin_p("phitop",phitop) if (is_master) write(*,*)trim(rname)// " phitop = ",phitop if (is_master) write(*,*)trim(rname)// "call molecular viscosity ?" callmolvis=.false. ! default value call getin_p("callmolvis",callmolvis) if (is_master) write(*,*)trim(rname)// " callmolvis = ",callmolvis if (is_master) write(*,*)trim(rname)// "call molecular diffusion ?" callmoldiff=.false. ! default value call getin_p("callmoldiff",callmoldiff) if (is_master) write(*,*)trim(rname)// " callmoldiff = ",callmoldiff ! 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)//": 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)//": 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 of radiatively active aerosols:" naerkind=0 ! default value call getin_p("naerkind",naerkind) if (is_master) write(*,*)trim(rname)//": naerkind = ",naerkind !*************************************************************** !! TRACERS options if (is_master)write(*,*)trim(rname)//& "call N2 condensation ?" n2cond=.true. ! default value call getin_p("n2cond",n2cond) if (is_master)write(*,*)trim(rname)//& " n2cond = ",n2cond if (is_master)write(*,*)trim(rname)//& "call N2 cloud condensation ?" condensn2=.false. ! default value call getin_p("condensn2",condensn2) if (is_master)write(*,*)trim(rname)//& "condensn2 = ",condensn2 if (is_master)write(*,*)trim(rname)//& "call no N2 frost formation?" no_n2frost=.false. ! default value call getin_p("no_n2frost",no_n2frost) if (is_master)write(*,*)trim(rname)//& "no_n2frost = ",no_n2frost if (is_master)write(*,*)trim(rname)//& "N2 condensation subtimestep?" nbsub=20 ! default value call getin_p("nbsub",nbsub) if (is_master)write(*,*)trim(rname)//& " nbsub = ",nbsub if (is_master)write(*,*)trim(rname)//& "Gravitationnal sedimentation ?" sedimentation=.true. ! default value call getin_p("sedimentation",sedimentation) if (is_master)write(*,*)trim(rname)//& " sedimentation = ",sedimentation if (is_master)write(*,*)trim(rname)//& "Compute glacial flow ?" glaflow=.false. ! default value call getin_p("glaflow",glaflow) if (is_master)write(*,*)trim(rname)//& " glaflow = ",glaflow if (is_master)write(*,*)trim(rname)//& "Compute methane cycle ?" methane=.false. ! default value call getin_p("methane",methane) if (is_master)write(*,*)trim(rname)//& " methane = ",methane condmetsurf=.true. ! default value call getin_p("condmetsurf",condmetsurf) if (is_master)write(*,*)trim(rname)//& " condmetsurf = ",condmetsurf if (is_master)write(*,*)trim(rname)//& "Compute methane clouds ?" metcloud=.false. ! default value call getin_p("metcloud",metcloud) if (is_master)write(*,*)trim(rname)//& " metcloud = ",metcloud if (is_master)write(*,*)trim(rname)//& "Compute CO cycle ?" carbox=.false. ! default value call getin_p("carbox",carbox) if (is_master)write(*,*)trim(rname)//& " carbox = ",carbox condcosurf=.true. ! default value call getin_p("condcosurf",condcosurf) if (is_master)write(*,*)trim(rname)//& " condcosurf = ",condcosurf if (is_master)write(*,*)trim(rname)//& "Compute CO clouds ?" monoxcloud=.false. ! default value call getin_p("monoxcloud",monoxcloud) if (is_master)write(*,*)trim(rname)//& " monoxcloud = ",monoxcloud if (is_master)write(*,*)trim(rname)//& "atmospheric redistribution (s):" tau_n2=1. ! default value call getin_p("tau_n2",tau_n2) if (is_master)write(*,*)trim(rname)//& " tau_n2 = ",tau_n2 tau_ch4=1.E7 ! default value call getin_p("tau_ch4",tau_ch4) if (is_master)write(*,*)trim(rname)//& " tau_ch4 = ",tau_ch4 tau_co=1. ! default value call getin_p("tau_co",tau_co) if (is_master)write(*,*)trim(rname)//& " tau_co = ",tau_co tau_prechaze=1. ! default value call getin_p("tau_prechaze",tau_prechaze) if (is_master)write(*,*)trim(rname)//& " tau_prechaze = ",tau_prechaze if (is_master)write(*,*)trim(rname)//& "Day fraction for limited cold trap in SP?" dayfrac=0. ! default value call getin_p("dayfrac",dayfrac) if (is_master)write(*,*)trim(rname)//& " dayfrac = ",dayfrac thresh_non2=0. ! default value call getin_p("thresh_non2",thresh_non2) if (is_master)write(*,*)trim(rname)//& " thresh_non2 = ",thresh_non2 ! use Pluto.old routines if (is_master) write(*,*) trim(rname)//": use vdifc from old Pluto ?" oldplutovdifc=.false. ! default value call getin_p("oldplutovdifc",oldplutovdifc) if (is_master) write(*,*) trim(rname)//": oldplutovdifc = ",oldplutovdifc if (is_master) write(*,*) trim(rname)//& ": call pluto.old correlated-k radiative transfer ?" oldplutocorrk=.false. ! default value call getin_p("oldplutocorrk",oldplutocorrk) if (is_master) write(*,*) trim(rname)//": oldplutocorrk = ",oldplutocorrk if (is_master) write(*,*) trim(rname)//& ": call pluto.old sedimentation ?" oldplutosedim=.false. ! default value call getin_p("oldplutosedim",oldplutosedim) if (is_master) write(*,*) trim(rname)//": oldplutosedim = ",oldplutosedim !*************************************************************** !! Haze options 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 hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed call getin_p("hazeconservch4",hazeconservch4) if (is_master)write(*,*)trim(rname)//& "hazconservch4 = ",hazeconservch4 if (is_master)write(*,*)trim(rname)//& "Production of haze (fast version) ?" fasthaze=.false. ! default value call getin_p("fasthaze",fasthaze) if (is_master)write(*,*)trim(rname)//& "fasthaze = ",fasthaze if (is_master)write(*,*)trim(rname)//& "Add source of haze ?" source_haze=.false. ! default value call getin_p("source_haze",source_haze) if (is_master)write(*,*)trim(rname)//& " source_haze = ",source_haze mode_hs=0 ! mode haze source default value call getin_p("mode_hs",mode_hs) if (is_master)write(*,*)trim(rname)//& " mode_hs",mode_hs kfix=1 ! default value call getin_p("kfix",kfix) if (is_master)write(*,*)trim(rname)//& "levels kfix",kfix fracsource=0.1 ! default value call getin_p("fracsource",fracsource) if (is_master)write(*,*)trim(rname)//& " fracsource",fracsource latsource=30. ! default value call getin_p("latsource",latsource) if (is_master)write(*,*)trim(rname)//& " latsource",latsource lonsource=180. ! default value call getin_p("lonsource",lonsource) if (is_master)write(*,*)trim(rname)//& " lonsource",lonsource if (is_master)write(*,*)trim(rname)//& "Radiatively active haze ?" aerohaze=.false. ! default value call getin_p("aerohaze",aerohaze) if (is_master)write(*,*)trim(rname)//& "aerohaze = ",aerohaze if (is_master)write(*,*)trim(rname)//& "Haze monomer radius ?" rad_haze=10.e-9 ! default value call getin_p("rad_haze",rad_haze) if (is_master)write(*,*)trim(rname)//& "rad_haze = ",rad_haze if (is_master)write(*,*)trim(rname)//& "fractal particle ?" fractal=.false. ! default value call getin_p("fractal",fractal) if (is_master)write(*,*)trim(rname)//& "fractal = ",fractal nb_monomer=10 ! default value call getin_p("nb_monomer",nb_monomer) if (is_master)write(*,*)trim(rname)//& "nb_monomer = ",nb_monomer if (is_master)write(*,*)trim(rname)//& "fixed gaseous CH4 mixing ratio for rad transfer ?" ch4fix=.false. ! default value call getin_p("ch4fix",ch4fix) if (is_master)write(*,*)trim(rname)//& " ch4fix = ",ch4fix if (is_master)write(*,*)trim(rname)//& "fixed gaseous CH4 mixing ratio for rad transfer ?" vmrch4fix=0.5 ! default value call getin_p("vmrch4fix",vmrch4fix) if (is_master)write(*,*)trim(rname)//& " vmrch4fix = ",vmrch4fix vmrch4_proffix=.false. ! default value call getin_p("vmrch4_proffix",vmrch4_proffix) if (is_master)write(*,*)trim(rname)//& " vmrch4_proffix = ",vmrch4_proffix if (is_master)write(*,*)trim(rname)//& "call specific cooling for radiative transfer ?" cooling=.false. ! default value call getin_p("cooling",cooling) if (is_master)write(*,*)trim(rname)//& " cooling = ",cooling if (is_master)write(*,*)trim(rname)//& "NLTE correction for methane heating rates?" nlte=.false. ! default value call getin_p("nlte",nlte) if (is_master)write(*,*)trim(rname)//& " nlte = ",nlte strobel=.false. ! default value call getin_p("strobel",strobel) if (is_master)write(*,*)trim(rname)//& " strobel = ",strobel if (is_master)write(*,*)trim(rname)//& "fixed radius profile from txt file ?" haze_radproffix=.false. ! default value call getin_p("haze_radproffix",haze_radproffix) if (is_master)write(*,*)trim(rname)//& "haze_radproffix = ",haze_radproffix if (is_master)write(*,*)trim(rname)//& "fixed MMR profile from txt file ?" haze_proffix=.false. ! default value call getin_p("haze_proffix",haze_proffix) if (is_master)write(*,*)trim(rname)//& "haze_proffix = ",haze_proffix if (is_master)write(*,*)trim(rname)//& "Number mix ratio of haze particles for co clouds:" Nmix_co=100000. ! default value call getin_p("Nmix_co",Nmix_co) if (is_master)write(*,*)trim(rname)//& " Nmix_co = ",Nmix_co if (is_master)write(*,*)trim(rname)//& "Number mix ratio of haze particles for ch4 clouds:" Nmix_ch4=100000. ! default value call getin_p("Nmix_ch4",Nmix_ch4) if (is_master)write(*,*)trim(rname)//& " Nmix_ch4 = ",Nmix_ch4 !*************************************************************** !! Surface properties !*********** N2 ********************************* if (is_master)write(*,*)trim(rname)//& "Mode for changing N2 albedo" mode_n2=0 ! default value call getin_p("mode_n2",mode_n2) if (is_master)write(*,*)trim(rname)//& " mode_n2 = ",mode_n2 thres_n2ice=1. ! default value call getin_p("thres_n2ice",thres_n2ice) if (is_master)write(*,*)trim(rname)//& " thres_n2ice = ",thres_n2ice if (is_master)write(*,*)trim(rname)//& "Diff of N2 albedo with thickness" deltab=0. ! default value call getin_p("deltab",deltab) if (is_master)write(*,*)trim(rname)//& " deltab = ",deltab if (is_master)write(*,*)trim(rname)//& "albedo N2 beta " alb_n2b=0.7 ! default value call getin_p("alb_n2b",alb_n2b) if (is_master)write(*,*)trim(rname)//& " alb_n2b = ",alb_n2b if (is_master)write(*,*)trim(rname)//& "albedo N2 alpha " alb_n2a=0.7 ! default value call getin_p("alb_n2a",alb_n2a) if (is_master)write(*,*)trim(rname)//& " alb_n2a = ",alb_n2a if (is_master)write(*,*)trim(rname)//& "emis N2 beta " emis_n2b=0.7 ! default value call getin_p("emis_n2b",emis_n2b) if (is_master)write(*,*)trim(rname)//& " emis_n2b = ",emis_n2b if (is_master)write(*,*)trim(rname)//& "emis N2 alpha " emis_n2a=0.7 ! default value call getin_p("emis_n2a",emis_n2a) if (is_master)write(*,*)trim(rname)//& " emis_n2a = ",emis_n2a !*********** CH4 ********************************* if (is_master)write(*,*)trim(rname)//& "Mode for changing CH4 albedo" mode_ch4=0 ! default value call getin_p("mode_ch4",mode_ch4) if (is_master)write(*,*)trim(rname)//& " mode_ch4 = ",mode_ch4 feedback_met=0 ! default value call getin_p("feedback_met",feedback_met) if (is_master)write(*,*)trim(rname)//& " feedback_met = ",feedback_met thres_ch4ice=1. ! default value, mm call getin_p("thres_ch4ice",thres_ch4ice) if (is_master)write(*,*)trim(rname)//& " thres_ch4ice = ",thres_ch4ice fdch4_latn=-1. fdch4_lats=0. fdch4_lone=0. fdch4_lonw=-1. fdch4_depalb=0.5 fdch4_finalb=0.95 fdch4_maxalb=0.99 fdch4_ampl=200. fdch4_maxice=100. call getin_p("fdch4_latn",fdch4_latn) call getin_p("fdch4_lats",fdch4_lats) call getin_p("fdch4_lone",fdch4_lone) call getin_p("fdch4_lonw",fdch4_lonw) call getin_p("fdch4_depalb",fdch4_depalb) call getin_p("fdch4_finalb",fdch4_finalb) call getin_p("fdch4_maxalb",fdch4_maxalb) call getin_p("fdch4_maxice",fdch4_maxice) call getin_p("fdch4_ampl",fdch4_ampl) if (is_master)write(*,*)trim(rname)//& " Values for albedo feedback = ",fdch4_latn,& fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,& fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl if (is_master)write(*,*)trim(rname)//& "Latitude for diff albedo methane" metlateq=25. ! default value call getin_p("metlateq",metlateq) if (is_master)write(*,*)trim(rname)//& " metlateq = ",metlateq if (is_master)write(*,*)trim(rname)//& "Ls1 and Ls2 for change of ch4 albedo" metls1=-1. ! default value metls2=-2. ! default value call getin_p("metls1",metls1) call getin_p("metls2",metls2) if (is_master)write(*,*)trim(rname)//& "albedo CH4 " alb_ch4=0.5 ! default value call getin_p("alb_ch4",alb_ch4) if (is_master)write(*,*)trim(rname)//& " alb_ch4 = ",alb_ch4 if (is_master)write(*,*)trim(rname)//& "albedo equatorial CH4 " alb_ch4_eq=alb_ch4 ! default value call getin_p("alb_ch4_eq",alb_ch4_eq) if (is_master)write(*,*)trim(rname)//& " alb_ch4_eq = ",alb_ch4_eq if (is_master)write(*,*)trim(rname)//& "albedo south hemis CH4 " alb_ch4_s=alb_ch4 ! default value call getin_p("alb_ch4_s",alb_ch4_s) if (is_master)write(*,*)trim(rname)//& " alb_ch4_s = ",alb_ch4_s if (is_master)write(*,*)trim(rname)//& "emis CH4 " emis_ch4=0.5 ! default value call getin_p("emis_ch4",emis_ch4) if (is_master)write(*,*)trim(rname)//& " emis_ch4 = ",emis_ch4 if (is_master)write(*,*)trim(rname)//& "CH4 lag for n2 sublimation limitation" ch4lag=.false. ! default value latlag=-90. ! default value vmrlag=1. ! default value call getin_p("ch4lag",ch4lag) call getin_p("latlag",latlag) call getin_p("vmrlag",vmrlag) if (is_master)write(*,*)trim(rname)//& " ch4lag = ",ch4lag if (is_master)write(*,*)trim(rname)//& " latlag = ",latlag if (is_master)write(*,*)trim(rname)//& " vmrlag = ",vmrlag if (is_master)write(*,*)trim(rname)//& "Max temperature for surface ?" tsurfmax=.false. ! default value albmin_ch4=0.3 ! default value call getin_p("tsurfmax",tsurfmax) call getin_p("albmin_ch4",albmin_ch4) if (is_master)write(*,*)trim(rname)//& " tsurfmax = ",tsurfmax if (is_master)write(*,*)trim(rname)//& " albmin_ch4 = ",albmin_ch4 if (is_master)write(*,*)trim(rname)//& "fixed gaseous CH4 mixing ratio for rad transfer ?" ch4fix=.false. ! default value call getin_p("ch4fix",ch4fix) if (is_master)write(*,*)trim(rname)//& " ch4fix = ",ch4fix if (is_master)write(*,*)trim(rname)//& "fixed gaseous CH4 mixing ratio for rad transfer ?" vmrch4fix=0.5 ! default value call getin_p("vmrch4fix",vmrch4fix) if (is_master)write(*,*)trim(rname)//& " vmrch4fix = ",vmrch4fix vmrch4_proffix=.false. ! default value call getin_p("vmrch4_proffix",vmrch4_proffix) if (is_master)write(*,*)trim(rname)//& " vmrch4_proffix = ",vmrch4_proffix !*********** CO ********************************* if (is_master)write(*,*)trim(rname)//& "albedo CO " alb_co=0.4 ! default value call getin_p("alb_co",alb_co) if (is_master)write(*,*)trim(rname)//& " alb_co = ",alb_co thres_coice=1. ! default value, mm call getin_p("thres_coice",thres_coice) if (is_master)write(*,*)trim(rname)//& " thres_coice = ",thres_coice if (is_master)write(*,*)trim(rname)//& "emis CO " emis_co=0.4 ! default value call getin_p("emis_co",emis_co) if (is_master)write(*,*)trim(rname)//& " emis_co = ",emis_co !*********** THOLINS ********************************* if (is_master)write(*,*)trim(rname)//& "Mode for changing tholins albedo/emis" mode_tholins=0 ! default value call getin_p("mode_tholins",mode_tholins) if (is_master)write(*,*)trim(rname)//& " mode_tholins = ",mode_tholins if (is_master)write(*,*)trim(rname)//& "albedo tho " alb_tho=0.1 ! default value call getin_p("alb_tho",alb_tho) if (is_master)write(*,*)trim(rname)//& " alb_tho = ",alb_tho if (is_master)write(*,*)trim(rname)//& "albedo tho eq" alb_tho_eq=0.1 ! default value call getin_p("alb_tho_eq",alb_tho_eq) if (is_master)write(*,*)trim(rname)//& " alb_tho_eq = ",alb_tho_eq if (is_master)write(*,*)trim(rname)//& "emis tholins " emis_tho=1. ! default value call getin_p("emis_tho",emis_tho) if (is_master)write(*,*)trim(rname)//& " emis_tho = ",emis_tho if (is_master)write(*,*)trim(rname)//& "emis tholins eq" emis_tho_eq=1. ! default value call getin_p("emis_tho_eq",emis_tho_eq) if (is_master)write(*,*)trim(rname)//& " emis_tho_eq = ",emis_tho_eq if (is_master)write(*,*)trim(rname)//& "Latitude for diff albedo tholins" tholateq=25. ! default value call getin_p("tholateq",tholateq) if (is_master)write(*,*)trim(rname)//& " tholateq = ",tholateq tholatn=-1. tholats=0. tholone=0. tholonw=-1. alb_tho_spe=0.1 ! default value emis_tho_spe=1. ! default value call getin_p(" tholatn",tholatn) call getin_p(" tholats",tholats) call getin_p(" tholone",tholone) call getin_p(" tholonw",tholonw) if (is_master)write(*,*)trim(rname)//& " Values for special tholins albedo = ",tholatn,& tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe if (is_master)write(*,*)trim(rname)//& "Specific albedo" spelon1=-180. ! default value spelon2=180. ! default value spelat1=-90. ! default value spelat2=90. ! default value specalb=.false. ! default value if (is_master)write(*,*)trim(rname)//& "albedo/emis spe " albspe=0.1 ! default value emispe=1. ! default value call getin_p("spelon1",spelon1) call getin_p("spelon2",spelon2) call getin_p("spelat1",spelat1) call getin_p("spelat2",spelat2) call getin_p("specalb",specalb) call getin_p("albspe",albspe) call getin_p("emispe",emispe) if (is_master)write(*,*)trim(rname)//& " specific = ",specalb !********************** TI ***************************** if (is_master)write(*,*)trim(rname)//& "Change TI with time" changeti=.false. ! default value call getin_p("changeti",changeti) if (is_master)write(*,*)trim(rname)//& " changeti = ",changeti changetid=.false. ! default value for diurnal TI call getin_p("changetid",changetid) if (is_master)write(*,*)trim(rname)//& " changetid = ",changetid if (is_master)write(*,*)trim(rname)//& "IT N2 " ITN2=800. ! default value call getin_p("ITN2",ITN2) if (is_master)write(*,*)trim(rname)//& " ITN2 = ",ITN2 ITN2d=20. ! default value call getin_p("ITN2d",ITN2d) if (is_master)write(*,*)trim(rname)//& " ITN2d = ",ITN2d if (is_master)write(*,*)trim(rname)//& "IT CH4" ITCH4=800. ! default value call getin_p("ITCH4",ITCH4) if (is_master)write(*,*)trim(rname)//& " ITCH4 = ",ITCH4 ITCH4d=20. ! default value call getin_p("ITCH4d",ITCH4d) if (is_master)write(*,*)trim(rname)//& " ITCH4d = ",ITCH4d if (is_master)write(*,*)trim(rname)//& "IT H2O" ITH2O=800. ! default value call getin_p("ITH2O",ITH2O) if (is_master)write(*,*)trim(rname)//& " ITH2O = ",ITH2O ITH2Od=20. ! default value call getin_p("ITH2Od",ITH2Od) if (is_master)write(*,*)trim(rname)//& " ITH2Od = ",ITH2Od !************** COOLING *************** alpha_top=5e-11 ! default value call getin_p("alpha_top",alpha_top) if (is_master)write(*,*)trim(rname)//& " alpha_top = ",alpha_top pref=0.02 ! default value call getin_p("pref",pref) if (is_master)write(*,*)trim(rname)//& " pref = ",pref deltap=0.1 ! default value call getin_p("deltap",deltap) if (is_master)write(*,*)trim(rname)//& " deltap = ",deltap !================================= ! 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)//& ": 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)//": Generic Condensation of tracers ?" generic_condensation=.false. !default value call getin_p("generic_condensation",generic_condensation) if (is_master) write(*,*)trim(rname)//": generic_condensation = ",generic_condensation 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)//": N2 ice albedo ?" albedon2ice=0.5 ! default value call getin_p("albedon2ice",albedon2ice) if (is_master) write(*,*)trim(rname)//": albedon2ice = ",albedon2ice 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)//": 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 ! Test of incompatibility: if ((.not.tracer).and.(haze)) then call abort_physic(rname, 'if haze are on, tracers must be on!', 1) endif if (haze_proffix.and.sedimentation) then call abort_physic(rname, 'if haze profile is set, sedimentation must be deactivated', 1) endif if (callmolvis.and..not.callconduct) then call abort_physic(rname, 'if callmolvis is set, callconduct must be true', 1) endif if (glaflow.and..not.fast) then call abort_physic(rname, 'if glaflow is set, fast must be true', 1) endif if (paleo.and..not.fast) then call abort_physic(rname, 'if paleo is set, fast must be true', 1) endif if ((haze_proffix.or.haze_radproffix).and..not.aerohaze) then call abort_physic(rname, 'for now, haze/rad proffix only works w aerohaze=T', 1) endif if (condcosurf.and.no_n2frost) then call abort_physic(rname, "CO surface condensation and no_n2frost are both active which may not be relevant", 1) end if 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) ! initialize variables and allocate arrays in radcommon_h call ini_radcommon_h(naerkind) ! allocate "comsoil_h" arrays call ini_comsoil_h(ngrid) END SUBROUTINE inifis END MODULE inifis_mod