SUBROUTINE meso_inifis(ngrid,nlayer,nq, $ wdt, $ wday_ini,wdaysec, $ wappel_phys,wecri_phys, $ plat,plon,parea, $ prad,pg,pr,pcpp, $ womeg,wmugaz, $ wyear_day,wperiheli,waphelie,wperi_day,wobliquit, $ wz0,wemin_turb,wlmixmin, $ wemissiv,wemissiceN,wemissiceS,walbediceN,walbediceS, $ wiceradiusN,wiceradiusS,wdtemisiceN,wdtemisiceS, $ walbedodat,wphisfi,wvolcapa, $ wzmea,wzstd,wzsig,wzgam,wzthe, $ wtheta,wpsi) c======================================================================= c c CAREFUL: THIS IS A VERSION TO BE USED WITH WRF !!! c c ... CHECK THE ****WRF lines c c======================================================================= ! !======================================================================= ! ! 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 ! adapted to the WRF use - Aymeric Spiga - Jan 2009 - Jan 2007 ! ! ! 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' !! **WRF a new stuff to be checked 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" #include "dimradmars.h" !!new #include "yomaer.h" !!new -- prob. pour remplir les proprietes #include "datafile.h" #include "meso_slope.h" #include "comsoil.h" !!**WRF -- needed to fill volcapa REAL prad,pg,pr,pcpp,pdaysec INTEGER ngrid,nlayer,nq REAL womeg,wmugaz,wdaysec REAL wyear_day,wperiheli,waphelie,wperi_day,wobliquit REAL wz0,wemin_turb,wlmixmin REAL wemissiv,wemissiceN,wemissiceS,walbediceN,walbediceS REAL wiceradiusN,wiceradiusS,wdtemisiceN,wdtemisiceS REAL walbedodat(ngrid),wphisfi(ngrid) REAL wzmea(ngrid),wzstd(ngrid),wzsig(ngrid) REAL wzgam(ngrid),wzthe(ngrid) REAL wtheta(ngrid),wpsi(ngrid) REAL wvolcapa REAL wdt REAL plat(ngrid),plon(ngrid),parea(ngridmx) integer wday_ini INTEGER ig,ierr INTEGER wecri_phys, wappel_phys ! EXTERNAL iniorbit,orbite EXTERNAL SSUM REAL SSUM CHARACTER ch1*12 CHARACTER ch80*80 #ifdef MESOSCALE ! logical chem, h2o ! chem = .false. ! h2o = .false. c ****WRF c c------------------------------------------------------ c Fill some parameters in the 'include' files c >> Do part of the job previously done by phyetat0.F c >> Complete list of parameters is found in tabfi.F c------------------------------------------------------ c c Values are defined in the module_model_constants.F WRF routine c ! in 'comcstfi.h' rad=prad cpp=pcpp g=pg r=pr rcp=r/cpp daysec=wdaysec omeg=womeg mugaz=wmugaz print*,"check: rad,cpp,g,r,rcp,daysec,omeg,mugaz" print*,rad,cpp,g,r,rcp,daysec,omeg,mugaz ! in 'planet.h' year_day=wyear_day periheli=wperiheli aphelie=waphelie peri_day=wperi_day obliquit=wobliquit z0=wz0 emin_turb=wemin_turb lmixmin=wlmixmin print*,"check: year_day,periheli,aphelie,peri_day,obliquit" print*,year_day,periheli,aphelie,peri_day,obliquit print*,"check: z0,emin_turb,lmixmin" print*,z0,emin_turb,lmixmin ! in 'surfdat.h' emissiv=wemissiv emisice(1)=wemissiceN emisice(2)=wemissiceS albedice(1)=walbediceN albedice(2)=walbediceS iceradius(1)=wiceradiusN iceradius(2)=wiceradiusS dtemisice(1)=wdtemisiceN dtemisice(2)=wdtemisiceS print*,"check: emissiv,emisice,albedice,iceradius,dtemisice" print*,emissiv,emisice,albedice,iceradius,dtemisice c c Values are defined in the WPS processing c albedodat(:)=walbedodat(:) !!!!! ***WRF inertiedat was moved, new physics !! !inertiedat(:)=winertiedat(:) phisfi(:)=wphisfi(:) print*,"check: albedodat(1),phisfi(1)" print*,albedodat(1),phisfi(1) print*,"check: albedodat(end),phisfi(end)" print*,albedodat(ngrid),phisfi(ngrid) ! NB: usually, gravity wave scheme is useless in mesoscale modeling ! NB: we however keep the option for coarse grid case ... zmea(:)=wzmea(:) zstd(:)=wzstd(:) zsig(:)=wzsig(:) zgam(:)=wzgam(:) zthe(:)=wzthe(:) print*,"check: gw param" print*,zmea(1),zmea(ngrid) print*,zstd(1),zstd(ngrid) print*,zsig(1),zsig(ngrid) print*,zgam(1),zgam(ngrid) print*,zthe(1),zthe(ngrid) ! ! in meso_slope.h ! theta_sl(:)=wtheta(:) psi_sl(:)=wpsi(:) print*,"check: theta_sl(1),psi_sl(1)" print*,theta_sl(1),psi_sl(1) print*,"check: theta_sl(end),psi_sl(end)" print*,theta_sl(ngrid),psi_sl(ngrid) ! ! in comsoil.h ! volcapa=wvolcapa print*,"check: volcapa" print*,volcapa c ****WRF ! rad=prad ! daysec=pdaysec ! dtphys=ptimestep ! cpp=pcpp ! g=pg ! r=pr ! rcp=r/cpp ! -------------------------------------------------------- ! 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: Parameters for the physics (callphys.def)' PRINT*,'--------------------------------------------' write(*,*) "Directory where external input files are:" !datafile="/u/forget/WWW/datagcm/datafile" ! NB: default 'datafile' is set in datafile.h call getin("datadir",datafile) write(*,*) " datafile = ",trim(datafile) 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(*,*) "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 callstats=.false. ! default value call getin("callstats",callstats) write(*,*) " callstats = ",callstats write(*,*) "Save EOF profiles in file 'profiles' for ", & "Climate Database?" calleofdump=.false. ! default value call getin("calleofdump",calleofdump) write(*,*) " calleofdump = ",calleofdump write(*,*) "Dust scenario: 1=constant dust (read from startfi", & " or set as tauvis); 2=Viking scenario; =3 MGS scenario,", & "=4 Mars Year 24 from TES assimilation, ", & "=24,25 or 26 :Mars Year 24,25 or 26 from TES assimilation" iaervar=3 ! default value call getin("iaervar",iaervar) write(*,*) " iaervar = ",iaervar write(*,*) "Reference (visible) dust opacity at 700 Pa ", & "(matters only if iaervar=1)" ! NB: default value of tauvis is set/read in startfi.nc file call getin("tauvis",tauvis) write(*,*) " tauvis = ",tauvis write(*,*) "Dust vertical distribution:" write(*,*) "(=1 top set by topdustref parameter;", & " =2 Viking scenario; =3 MGS scenario)" 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 NLTE radiative schemes ?", & "(matters only if callrad=T)" callnlte=.false. ! default value call getin("callnlte",callnlte) write(*,*) " callnlte = ",callnlte write(*,*) "call CO2 NIR absorption ?", & "(matters only if callrad=T)" callnirco2=.false. ! default value call getin("callnirco2",callnirco2) write(*,*) " callnirco2 = ",callnirco2 write(*,*) "call turbulent vertical diffusion ?" calldifv=.true. ! default value call getin("calldifv",calldifv) write(*,*) " calldifv = ",calldifv write(*,*) "call thermals ?" calltherm=.false. ! default value call getin("calltherm",calltherm) write(*,*) " calltherm = ",calltherm write(*,*) "output thermal diagnostics ?" outptherm=.false. ! default value call getin("outptherm",outptherm) write(*,*) " outptherm = ",outptherm write(*,*) "call convective adjustment ?" calladj=.true. ! default value call getin("calladj",calladj) write(*,*) " calladj = ",calladj if (calltherm .and. (.not. calladj)) then print*,'Convadj has to be activated when using thermals' stop endif write(*,*) "call CO2 condensation ?" callcond=.true. ! default value call getin("callcond",callcond) write(*,*) " callcond = ",callcond 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(*,*)"Output of the exchange coefficient mattrix ?", & "(for diagnostics only)" callg2d=.false. ! default value call getin("callg2d",callg2d) write(*,*)" callg2d = ",callg2d write(*,*)"Rayleigh scattering : (should be .false. for now)" rayleigh=.false. call getin("rayleigh",rayleigh) write(*,*)" rayleigh = ",rayleigh ! TRACERS: write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)" dustbin=0 ! default value call getin("dustbin",dustbin) write(*,*)" dustbin = ",dustbin write(*,*)"Radiatively active dust ? (matters if dustbin>0)" active=.false. ! default value call getin("active",active) write(*,*)" active = ",active ! 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 submicron=.false. ! default value call getin("submicron",submicron) write(*,*)" submicron = ",submicron ! Test of incompatibility: ! if doubleq is used, then dustbin should be 2 if (doubleq.and.(dustbin.ne.2)) then print*,'if doubleq is used, then dustbin should be 2' stop endif if (doubleq.and.submicron.and.(nqmx.LT.3)) then print*,'If doubleq is used with a submicron tracer,' print*,' then the number of tracers has to be' print*,' larger than 3.' 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(*,*) "Radiatively active transported atmospheric ", & "water ice ?" activice=.false. ! default value call getin("activice",activice) write(*,*) " activice = ",activice write(*,*) "Compute water cycle ?" water=.false. ! default value call getin("water",water) write(*,*) " water = ",water ! Test of incompatibility: if (activice.and..not.water) then print*,'if activice is used, water should be used too' stop endif if (water.and..not.tracer) then print*,'if water is used, tracer should be used too' stop endif ! Test of incompatibility: write(*,*) "Permanent water caps at poles ?", & " .true. is RECOMMENDED" write(*,*) "(with .true., North cap is a source of water ", & "and South pole is a cold trap)" caps=.true. ! default value call getin("caps",caps) write(*,*) " caps = ",caps 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) c ****WRF c***************************************************** c Since it comes from WRF settings, we have to c fill dtphys in the include file c It must be set now, because it is used afterwards c c Opportunity is taken to fill ecri_phys as well c***************************************************** dtphys=wdt*float(wappel_phys) print*,'Physical timestep (s) ',dtphys ecri_phys=wecri_phys print*,'Physical frequency for writing ',ecri_phys c c ****WRF 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=1 ilwn=1 !2 ilwb=1 !2 linear=.true. ncouche=3 alphan=0.4 semi=0 ! BUT people working hard on the LW may want to read them in 'radia.def' ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ OPEN(99,file='radia.def',status='old',form='formatted' . ,iostat=ierr) IF(ierr.EQ.0) THEN write(*,*) 'inifis: Reading radia.def !!!' READ(99,fmt='(a)') ch1 READ(99,*) callemis WRITE(*,8000) ch1,callemis READ(99,fmt='(a)') ch1 READ(99,*) iradia WRITE(*,8001) ch1,iradia READ(99,fmt='(a)') ch1 READ(99,*) ilwd WRITE(*,8001) ch1,ilwd READ(99,fmt='(a)') ch1 READ(99,*) ilwn WRITE(*,8001) ch1,ilwn READ(99,fmt='(a)') ch1 READ(99,*) linear WRITE(*,8000) ch1,linear READ(99,fmt='(a)') ch1 READ(99,*) ncouche WRITE(*,8001) ch1,ncouche READ(99,fmt='(a)') ch1 READ(99,*) alphan WRITE(*,*) ch1,alphan READ(99,fmt='(a)') ch1 READ(99,*) ilwb WRITE(*,8001) ch1,ilwb READ(99,fmt='(a)') ch1 READ(99,'(l1)') callg2d WRITE(*,8000) ch1,callg2d READ(99,fmt='(a)') ch1 READ(99,*) semi WRITE(*,*) ch1,semi end if CLOSE(99) !----------------------------------------------------------------------- ! Some more initialization: ! ------------------------ ! in 'comgeomfi.h' 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) ! in 'comdiurn.h' 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: ! ------------------------------- ! Ehouarn: removed; as these tests are now done in initracer.F ! 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) h2o = .true. ! !! TESTS ! ! print*,'inifis: TRACERS:' ! write(*,*) " chem=",chem," h2o=",h2o !! write(*,*) " doubleq=",doubleq !! write(*,*) " dustbin=",dustbin ! ! if ((doubleq).and.(h2o).and. ! $ (chem)) then ! print*,' 2 dust tracers (doubleq)' ! print*,' 1 water vapour tracer' ! print*,' 1 water ice tracer' ! print*,nqmx-4,' chemistry tracers' ! endif ! ! if ((doubleq).and.(h2o).and. ! $ .not.(chem)) then ! print*,' 2 dust tracers (doubleq)' ! print*,' 1 water vapour tracer' ! print*,' 1 water ice tracer' ! if (nqmx.LT.4) then ! print*,'nqmx should be at least equal to' ! print*,'4 with these options.' ! stop ! endif ! endif ! ! if (.not.(doubleq).and.(h2o).and. ! $ (chem)) 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. ! $ .not.(chem)) 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 ! write(*,*) "inifis: nqmx.lt.(dustbin+2)" ! stop ! endif ! endif ! ! endif ! of if (tracer) ! ! RETURN #endif END