SUBROUTINE inifis(ngrid,nlayer, $ day_ini,pdaysec,ptimestep, $ plat,plon,parea, $ prad,pg,pr,pcpp) ! !======================================================================= ! ! purpose: ! ------- ! ! Initialisation for the physical parametrisations of the LMD ! martian atmospheric general circulation modele. ! ! 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: ! ------------- ! to use 'getin' USE ioipsl_getincom IMPLICIT NONE #include "dimensions.h" #include "dimphys.h" #include "planete.h" #include "comcstfi.h" #include "comsaison.h" #include "comdiurn.h" #include "comgeomfi.h" #include "callkeys.h" #include "surfdat.h" REAL prad,pg,pr,pcpp,pdaysec,ptimestep INTEGER ngrid,nlayer REAL plat(ngrid),plon(ngrid),parea(ngridmx) integer day_ini INTEGER ig,ierr EXTERNAL iniorbit,orbite EXTERNAL SSUM REAL SSUM CHARACTER ch1*12 CHARACTER ch80*80 logical chem, h2o logical :: parameter, doubleq=.false. rad=prad daysec=pdaysec dtphys=ptimestep cpp=pcpp g=pg r=pr rcp=r/cpp avocado = 6.02214179e23 ! added by RW ! -------------------------------------------------------- ! The usual Tests ! -------------- IF (nlayer.NE.nlayermx) THEN PRINT*,'STOP in inifis' PRINT*,'Probleme de dimensions :' PRINT*,'nlayer = ',nlayer PRINT*,'nlayermx = ',nlayermx STOP ENDIF IF (ngrid.NE.ngridmx) THEN PRINT*,'STOP in inifis' PRINT*,'Probleme de dimensions :' PRINT*,'ngrid = ',ngrid PRINT*,'ngridmx = ',ngridmx STOP ENDIF ! -------------------------------------------------------------- ! Reading the "callphys.def" file controlling some key options ! -------------------------------------------------------------- ! check that 'callphys.def' file is around OPEN(99,file='callphys.def',status='old',form='formatted' & ,iostat=ierr) CLOSE(99) IF(ierr.EQ.0) THEN PRINT* PRINT* PRINT*,'--------------------------------------------' PRINT*,' inifis: Parametres pour la physique (callphys.def)' PRINT*,'--------------------------------------------' write(*,*) "Run with or without tracer transport ?" tracer=.false. ! default value call getin("tracer",tracer) write(*,*) " tracer = ",tracer write(*,*) "Diurnal cycle ?" write(*,*) "(if diurnal=false, diurnal averaged solar heating)" diurnal=.true. ! default value call getin("diurnal",diurnal) write(*,*) " diurnal = ",diurnal write(*,*) "Seasonal cycle ?" write(*,*) "(if season=false, Ls stays constant, to value ", & "set in 'start'" season=.true. ! default value call getin("season",season) write(*,*) " season = ",season write(*,*) "Tidally resonant rotation ?" tlocked=.false. ! default value call getin("tlocked",tlocked) write(*,*) "tlocked = ",tlocked ! Test of incompatibility: ! if tlocked, then diurnal should be false if (tlocked.and.diurnal) then print*,'If diurnal=true, we should turn off tlocked.' stop endif write(*,*) "Tidal resonance ratio ?" nres=0 ! default value call getin("nres",nres) write(*,*) "nres = ",nres write(*,*) "Write some extra output to the screen ?" lwrite=.false. ! default value call getin("lwrite",lwrite) write(*,*) " lwrite = ",lwrite write(*,*) "Save statistics in file stats.nc ?" callstats=.true. ! default value call getin("callstats",callstats) write(*,*) " callstats = ",callstats write(*,*) "Test energy conservation of model physics ?" enertest=.false. ! default value call getin("enertest",enertest) write(*,*) " enertest = ",enertest write(*,*) "Save EOF profiles in file 'profiles' for ", & "Climate Database?" calleofdump=.false. ! default value call getin("calleofdump",calleofdump) write(*,*) " calleofdump = ",calleofdump write(*,*) "Dust scenario:" iaervar=3 ! default value call getin("iaervar",iaervar) write(*,*) " iaervar = ",iaervar write(*,*) "Dust vertical distribution:" write(*,*) "(=1 Dust opt.deph read in startfi;", & " =2 Viking scenario; =3 MGS scenario,", & " =4 Mars Year 24 from TES assimilation)" iddist=3 ! default value call getin("iddist",iddist) write(*,*) " iddist = ",iddist write(*,*) "Dust top altitude (km). (Matters only if iddist=1)" topdustref= 90.0 ! default value call getin("topdustref",topdustref) write(*,*) " topdustref = ",topdustref write(*,*) "call radiative transfer ?" callrad=.true. ! default value call getin("callrad",callrad) write(*,*) " callrad = ",callrad write(*,*) "call correlated-k radiative transfer ?" corrk=.true. ! default value call getin("corrk",corrk) write(*,*) " corrk = ",corrk ! write(*,*) "call NLTE radiative schemes ?", ! & "(matters only if callrad=T)" ! callnlte=.false. ! default value ! call getin("callnlte",callnlte) ! write(*,*) " callnlte = ",callnlte write(*,*) "call gaseous absorption in the visible bands?", & "(matters only if callrad=T)" callgasvis=.false. ! default value call getin("callgasvis",callgasvis) write(*,*) " callgasvis = ",callgasvis write(*,*) "call turbulent vertical diffusion ?" calldifv=.true. ! default value call getin("calldifv",calldifv) write(*,*) " calldifv = ",calldifv write(*,*) "call convective adjustment ?" calladj=.true. ! default value call getin("calladj",calladj) write(*,*) " calladj = ",calladj write(*,*)"Gas is nonideal CO2 ?" nonideal=.false. call getin("nonideal",nonideal) write(*,*)" nonideal = ",nonideal ! Test of incompatibility if (enertest.and.nonideal) then print*,'Energy conservation calculations currently & assume ideal gas!' call abort endif write(*,*) "call CO2 condensation ?" co2cond=.true. ! default value call getin("co2cond",co2cond) write(*,*) " co2cond = ",co2cond write(*,*)"call thermal conduction in the soil ?" callsoil=.true. ! default value call getin("callsoil",callsoil) write(*,*) " callsoil = ",callsoil ! write(*,*)"call Lott's gravity wave/subgrid topography ", ! & "scheme ?" ! calllott=.true. ! default value ! call getin("calllott",calllott) ! write(*,*)" calllott = ",calllott write(*,*)"rad.transfer is computed every iradia", & " physical timestep" iradia=1 ! default value call getin("iradia",iradia) write(*,*)" iradia = ",iradia write(*,*)"Rayleigh scattering ?" rayleigh=.false. call getin("rayleigh",rayleigh) write(*,*)" rayleigh = ",rayleigh write(*,*)"Parametrized Earth-like ozone absorption ?" ozone=.false. call getin("ozone",ozone) write(*,*) " ozone = ",ozone write(*,*)"Output mean OLR in 1D?" meanOLR=.false. call getin("meanOLR",meanOLR) write(*,*)" meanOLR = ",meanOLR write(*,*)"Output spectral OLR in 3D?" specOLR=.false. call getin("specOLR",specOLR) write(*,*)" specOLR = ",specOLR write(*,*)"Default planetary temperature?" tplanet=215.0 call getin("tplanet",tplanet) write(*,*)" tplanet = ",tplanet write(*,*)"Which star?" startype=1 ! default value = Sol call getin("startype",startype) write(*,*)" startype = ",startype write(*,*)"Value of stellar flux at 1 AU?" Fat1AU=1356.0 ! default value = Sol today call getin("Fat1AU",Fat1AU) write(*,*)" Fat1AU = ",Fat1AU write(*,*)"Set temperature to 1 K above CO2 condensation?" nearco2cond=.false. call getin("nearco2cond",nearco2cond) write(*,*)" nearco2cond = ",nearco2cond ! TRACERS: write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)" dustbin=0 ! default value call getin("dustbin",dustbin) write(*,*)" dustbin = ",dustbin write(*,*)"Radiatively active aerosols?" aerofixed=.true. ! default value call getin("aerofixed",aerofixed) write(*,*)" aerofixed = ",aerofixed write(*,*)"Number mixing ratio of CO2 ice particles:" Nmix_co2=100000. ! default value call getin("Nmix_co2",Nmix_co2) write(*,*)" Nmix_co2 = ",Nmix_co2 write(*,*)"Number mixing ratio of H2O ice particles:" Nmix_h2o=10000000. ! default value call getin("Nmix_h2o",Nmix_h2o) write(*,*)" Nmix_h2o = ",Nmix_h2o write(*,*)"Is the variable gas species radiatively active?" varactive=.false. call getin("varactive",varactive) write(*,*)" varactive = ",varactive write(*,*)"Is the variable gas species distribution set?" varfixed=.false. call getin("varfixed",varfixed) write(*,*)" varfixed = ",varfixed write(*,*)"What is the saturation % of the variable species?" satval=0.8 call getin("satval",satval) write(*,*)" satval = ",satval ! write(*,*)"Radiatively active dust ? (matters if dustbin>0)" ! active=.false. ! default value ! call getin("active",active) ! write(*,*)" active = ",active ! Test of incompatibility: ! if no tracers, then aerofixed should be true if ((.not.tracer).and.(.not.aerofixed)) then print*,'if tracers are off, aerofixed must be ON!' stop endif ! Test of incompatibility: ! if varactive, then varfixed should be false if (varactive.and.varfixed) then print*,'if varactive, varfixed must be OFF!' stop endif ! Test of incompatibility: ! if active is used, then dustbin should be > 0 ! if (active.and.(dustbin.lt.1)) then ! print*,'if active is used, then dustbin should > 0' ! stop ! endif ! write(*,*)"use mass and number mixing ratios to predict", ! & " dust size ?" ! doubleq=.false. ! default value ! call getin("doubleq",doubleq) ! write(*,*)" doubleq = ",doubleq ! Test of incompatibility: ! if doubleq is used, then dustbin should be 1 ! if (doubleq.and.(dustbin.ne.1)) then ! print*,'if doubleq is used, then dustbin should be 1' ! stop ! endif ! write(*,*)"dust lifted by GCM surface winds ?" ! lifting=.false. ! default value ! call getin("lifting",lifting) ! write(*,*)" lifting = ",lifting ! Test of incompatibility: ! if lifting is used, then dustbin should be > 0 ! if (lifting.and.(dustbin.lt.1)) then ! print*,'if lifting is used, then dustbin should > 0' ! stop ! endif ! write(*,*)" dust lifted by dust devils ?" ! callddevil=.false. !default value ! call getin("callddevil",callddevil) ! write(*,*)" callddevil = ",callddevil ! Test of incompatibility: ! if dustdevil is used, then dustbin should be > 0 ! if (callddevil.and.(dustbin.lt.1)) then ! print*,'if dustdevil is used, then dustbin should > 0' ! stop ! endif ! write(*,*)"Dust scavenging by CO2 snowfall ?" ! scavenging=.false. ! default value ! call getin("scavenging",scavenging) ! write(*,*)" scavenging = ",scavenging ! Test of incompatibility: ! if scavenging is used, then dustbin should be > 0 ! if (scavenging.and.(dustbin.lt.1)) then ! print*,'if scavenging is used, then dustbin should > 0' ! stop ! endif write(*,*) "Gravitationnal sedimentation ?" sedimentation=.true. ! default value call getin("sedimentation",sedimentation) write(*,*) " sedimentation = ",sedimentation ! write(*,*) "includes water ice", ! & "(if true, 'water' must also be .true.)" ! iceparty=.false. ! default value ! call getin("iceparty",iceparty) ! write(*,*) " iceparty = ",iceparty ! write(*,*) "Radiatively active transported atmospheric ", ! & "water ice ?" ! activice=.false. ! default value ! call getin("activice",activice) ! write(*,*) " activice = ",activice ! Test of incompatibility: ! if activice is used, then iceparty should be used too ! if (activice.and..not.iceparty) then ! print*,'if activice is used, iceparty should be used too' ! stop ! endif write(*,*) "Compute water cycle ?" water=.false. ! default value call getin("water",water) write(*,*) " water = ",water write(*,*) "Include water condensation ?" watercond=.false. ! default value call getin("watercond",watercond) write(*,*) " watercond = ",watercond write(*,*) "Include water precipitation ?" waterrain=.false. ! default value call getin("waterrain",waterrain) write(*,*) " waterrain = ",waterrain write(*,*) "Precipitation threshold ?" rainthreshold=0.011 ! default value (Emmanuel 1997) call getin("rainthreshold",rainthreshold) write(*,*) " rainthreshold = ",rainthreshold ! Test of incompatibility: ! if watercond is used, then water should be used too if (watercond.and.(.not.watercond)) then print*,'if watercond is used, water should be used too' stop endif write(*,*) "photochemistry: include chemical species" photochem=.false. ! default value call getin("photochem",photochem) write(*,*) " photochem = ",photochem ! THERMOSPHERE write(*,*) "call thermosphere ?" callthermos=.false. ! default value call getin("callthermos",callthermos) write(*,*) " callthermos = ",callthermos write(*,*) " water included without cycle ", & "(only if water=.false.)" thermoswater=.false. ! default value call getin("thermoswater",thermoswater) write(*,*) " thermoswater = ",thermoswater write(*,*) "call thermal conduction ?", & " (only if callthermos=.true.)" callconduct=.false. ! default value call getin("callconduct",callconduct) write(*,*) " callconduct = ",callconduct write(*,*) "call EUV heating ?", & " (only if callthermos=.true.)" calleuv=.false. ! default value call getin("calleuv",calleuv) write(*,*) " calleuv = ",calleuv write(*,*) "call molecular viscosity ?", & " (only if callthermos=.true.)" callmolvis=.false. ! default value call getin("callmolvis",callmolvis) write(*,*) " callmolvis = ",callmolvis write(*,*) "call molecular diffusion ?", & " (only if callthermos=.true.)" callmoldiff=.false. ! default value call getin("callmoldiff",callmoldiff) write(*,*) " callmoldiff = ",callmoldiff write(*,*) "call thermospheric photochemistry ?", & " (only if callthermos=.true.)" thermochem=.false. ! default value call getin("thermochem",thermochem) write(*,*) " thermochem = ",thermochem write(*,*) "date for solar flux calculation:", & " (1985 < date < 2002)" write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)" solarcondate=1993.4 ! default value call getin("solarcondate",solarcondate) write(*,*) " solarcondate = ",solarcondate if (.not.callthermos) then if (thermoswater) then print*,'if thermoswater is set, callthermos must be true' stop endif if (callconduct) then print*,'if callconduct is set, callthermos must be true' stop endif if (calleuv) then print*,'if calleuv is set, callthermos must be true' stop endif if (callmolvis) then print*,'if callmolvis is set, callthermos must be true' stop endif if (callmoldiff) then print*,'if callmoldiff is set, callthermos must be true' stop endif if (thermochem) then print*,'if thermochem is set, callthermos must be true' stop endif endif ! Test of incompatibility: ! if photochem is used, then water should be used too if (photochem.and..not.water) then print*,'if photochem is used, water should be used too' stop endif ! if callthermos is used, then thermoswater should be used too ! (if water not used already) if (callthermos .and. .not.water) then if (callthermos .and. .not.thermoswater) then print*,'if callthermos is used, water or thermoswater & should be used too' stop endif endif PRINT*,'--------------------------------------------' PRINT* PRINT* ELSE write(*,*) write(*,*) 'Cannot read file callphys.def. Is it here ?' stop ENDIF 8000 FORMAT(t5,a12,l8) 8001 FORMAT(t5,a12,i8) 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* ! -------------------------------------------------------------- ! Managing the Longwave radiative transfer ! -------------------------------------------------------------- ! In most cases, the run just use the following values : ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ callemis=.true. ! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01 ilwd=10*int(daysec/dtphys) ilwn=2 linear=.true. ncouche=3 alphan=0.4 ilwb=2 semi=0 c$$$! BUT people working hard on the LW may want to read them in 'radia.def' c$$$! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c$$$ c$$$ OPEN(99,file='radia.def',status='old',form='formatted' c$$$ . ,iostat=ierr) c$$$ IF(ierr.EQ.0) THEN c$$$ write(*,*) 'inifis: Reading radia.def !!!' c$$$ READ(99,fmt='(a)') ch1 c$$$ READ(99,*) callemis c$$$ WRITE(*,8000) ch1,callemis c$$$ c$$$ READ(99,fmt='(a)') ch1 c$$$ READ(99,*) iradia c$$$ WRITE(*,8001) ch1,iradia c$$$ c$$$ READ(99,fmt='(a)') ch1 c$$$ READ(99,*) ilwd c$$$ WRITE(*,8001) ch1,ilwd c$$$ c$$$ READ(99,fmt='(a)') ch1 c$$$ READ(99,*) ilwn c$$$ WRITE(*,8001) ch1,ilwn c$$$ c$$$ READ(99,fmt='(a)') ch1 c$$$ READ(99,*) linear c$$$ WRITE(*,8000) ch1,linear c$$$ c$$$ READ(99,fmt='(a)') ch1 c$$$ READ(99,*) ncouche c$$$ WRITE(*,8001) ch1,ncouche c$$$ c$$$ READ(99,fmt='(a)') ch1 c$$$ READ(99,*) alphan c$$$ WRITE(*,*) ch1,alphan c$$$ c$$$ READ(99,fmt='(a)') ch1 c$$$ READ(99,*) ilwb c$$$ WRITE(*,8001) ch1,ilwb c$$$ c$$$ c$$$ READ(99,fmt='(a)') ch1 c$$$ READ(99,'(l1)') callg2d c$$$ WRITE(*,8000) ch1,callg2d c$$$ c$$$ READ(99,fmt='(a)') ch1 c$$$ READ(99,*) semi c$$$ WRITE(*,*) ch1,semi c$$$ end if c$$$ CLOSE(99) !----------------------------------------------------------------------- ! Some more initialization: ! ------------------------ CALL SCOPY(ngrid,plon,1,long,1) CALL SCOPY(ngrid,plat,1,lati,1) CALL SCOPY(ngrid,parea,1,area,1) totarea=SSUM(ngridmx,area,1) 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 pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h ! managing the tracers, and tests: ! ------------------------------- if(tracer) then ! when photochem is used, nqchem_min is the rank ! of the first chemical species ! Ehouarn: nqchem_min is now meaningless and no longer used ! nqchem_min = 1 if (photochem .or. callthermos) then chem = .true. end if if (water .or. thermoswater)then h2o = .true. else h2o = .false. endif ! TESTS print*,'inifis: TRACERS:' c$$$ if ((doubleq).and.(h2o).and. c$$$ $ (chem).and.(iceparty)) then c$$$ print*,' 2 dust tracers (doubleq)' c$$$ print*,' 1 water vapour tracer' c$$$ print*,' 1 water ice tracer' c$$$ print*,nqmx-4,' chemistry tracers' c$$$ endif c$$$ c$$$ if ((doubleq).and.(h2o).and. c$$$ $ (chem).and..not.(iceparty)) then c$$$ print*,' 2 dust tracers (doubleq)' c$$$ print*,' 1 water vapour tracer' c$$$ print*,nqmx-3,' chemistry tracers' c$$$ endif c$$$ c$$$ if ((doubleq).and.(h2o).and. c$$$ $ .not.(chem).and.(iceparty)) then c$$$ print*,' 2 dust tracers (doubleq)' c$$$ print*,' 1 water vapour tracer' c$$$ print*,' 1 water ice tracer' c$$$ if (nqmx.ne.4) then c$$$ print*,'nqmx should be 4 with these options.' c$$$ print*,'(or check callphys.def)' c$$$ stop c$$$ endif c$$$ endif c$$$ c$$$ if ((doubleq).and.(h2o).and. c$$$ $ .not.(chem).and..not.(iceparty)) then c$$$ print*,' 2 dust tracers (doubleq)' c$$$ print*,' 1 water vapour tracer' c$$$ if (nqmx.ne.3) then c$$$ print*,'nqmx should be 3 with these options...' c$$$ print*,'(or check callphys.def)' c$$$ stop c$$$ endif c$$$ endif c$$$ c$$$ if ((doubleq).and..not.(h2o)) then c$$$ print*,' 2 dust tracers (doubleq)' c$$$ if (nqmx.ne.2) then c$$$ print*,'nqmx should be 2 with these options...' c$$$ print*,'(or check callphys.def)' c$$$ stop c$$$ endif c$$$ endif if (.not.(doubleq).and.(h2o).and. ! $ (chem).and.(iceparty)) then $ (chem).and.(watercond)) then if (dustbin.gt.0) then print*,dustbin,' dust bins' endif print*,nqmx-2-dustbin,' chemistry tracers' print*,' 1 water vapour tracer' print*,' 1 water ice tracer' endif if (.not.(doubleq).and.(h2o).and. ! $ (chem).and..not.(iceparty)) then $ (chem).and..not.(watercond)) then if (dustbin.gt.0) then print*,dustbin,' dust bins' endif print*,nqmx-1-dustbin,' chemistry tracers' print*,' 1 water vapour tracer' endif if (.not.(doubleq).and.(h2o).and. ! $ .not.(chem).and.(iceparty)) then $ .not.(chem).and.(watercond)) then if (dustbin.gt.0) then print*,dustbin,' dust bins' endif print*,' 1 water vapour tracer' print*,' 1 water ice tracer' if (nqmx.gt.(dustbin+2)) then print*,'nqmx should be ',(dustbin+2), $ ' with these options...' print*,'(or check callphys.def)' endif if (nqmx.lt.(dustbin+2)) then print*,dustbin,' dust bins, but this should be ok I think.' ! stop endif endif if (.not.(doubleq).and.(h2o).and. ! $ .not.(chem).and..not.(iceparty)) then $ .not.(chem).and..not.(watercond)) then if (dustbin.gt.0) then print*,dustbin,' dust bins' endif print*,' 1 water vapour tracer' if (nqmx.gt.(dustbin+1)) then print*,'nqmx should be ',(dustbin+1), $ ' with these options...' print*,'(or check callphys.def)' if (nqmx.lt.(dustbin+1)) then stop endif endif endif ! if (.not.(doubleq).and..not.(h2o)) then ! if (dustbin.gt.0) then ! print*,dustbin,' dust bins' ! if (nqmx.ne.dustbin) then ! print*,'nqmx should be ',dustbin, ! $ ' with these options...' ! print*,'(or check callphys.def)' ! stop ! endif ! else print*,'dustbin=',dustbin, $ ': tracer should be F with these options...' $ ,'UNLESS you just want to move tracers around ' ! endif ! endif endif ! of if (tracer) RETURN END