source: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 @ 3533

Last change on this file since 3533 was 3497, checked in by debatzbr, 2 months ago

Add AC6H6 condensation in the microphysics
File size: 22.0 KB
RevLine 
[1524]1MODULE inifis_mod
2IMPLICIT NONE
[253]3
[1524]4CONTAINS
[253]5
[1524]6  SUBROUTINE inifis(ngrid,nlayer,nq, &
[1525]7             day_ini,pdaysec,nday,ptimestep, &
[1524]8             plat,plon,parea, &
9             prad,pg,pr,pcpp)
10
[1896]11  use init_print_control_mod, only: init_print_control
[1788]12  use radinc_h, only: ini_radinc_h
[1524]13  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
[1542]14  use comgeomfi_h, only: totarea, totarea_planet
[1538]15  use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
[1525]16  use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
17                              init_time, daysec, dtphys
[1524]18  use comcstfi_mod, only: rad, cpp, g, r, rcp, &
[1672]19                          mugaz, pi, avocado, kbol
[1524]20  use planete_mod, only: nres
21  use planetwide_mod, only: planetwide_sumval
22  use callkeys_mod
[1529]23  use mod_phys_lmdz_para, only : is_parallel
[1524]24
[135]25!=======================================================================
26!
27!   purpose:
28!   -------
29!
30!   Initialisation for the physical parametrisations of the LMD
[305]31!   Generic Model.
[135]32!
33!   author: Frederic Hourdin 15 / 10 /93
34!   -------
[3090]35!   modified: -Sebastien Lebonnois 11/06/2003 (new callphys.def)
36!             -Ehouarn Millour (oct. 2008) tracers are now identified
37!               by their names and may not be contiguously
38!               stored in the q(:,:,:,:) array
39!             -Ehouarn Millour (june 2009) use getin routine
40!               to load parameters
[3497]41!             -B. de Batz de Trenquelléon (2023) minor changes and
[3090]42!               new initialisations
[135]43!
44!
45!   arguments:
46!   ----------
47!
48!   input:
49!   ------
50!
51!    ngrid                 Size of the horizontal grid.
52!                          All internal loops are performed on that grid.
53!    nlayer                Number of vertical layers.
54!    pdayref               Day of reference for the simulation
55!    pday                  Number of days counted from the North. Spring
56!                          equinoxe.
57!
58!=======================================================================
59!
60!-----------------------------------------------------------------------
61!   declarations:
62!   -------------
[2366]63  use datafile_mod
[1524]64  use ioipsl_getin_p_mod, only: getin_p
65  IMPLICIT NONE
[1384]66
[135]67
68
[1524]69  REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
[1525]70  INTEGER,INTENT(IN) :: nday
[1524]71  INTEGER,INTENT(IN) :: ngrid,nlayer,nq
72  REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
73  integer,intent(in) :: day_ini
74  INTEGER ig,ierr
[135]75 
[1524]76  EXTERNAL iniorbit,orbite
77  EXTERNAL SSUM
78  REAL SSUM
[135]79 
[1896]80  ! Initialize flags lunout, prt_level, debug (in print_control_mod)
81  CALL init_print_control
82
[1524]83  ! initialize constants in comcstfi_mod
84  rad=prad
85  cpp=pcpp
86  g=pg
87  r=pr
88  rcp=r/cpp
[2078]89  mugaz=8.314*1000./pr ! dummy init
[1524]90  pi=2.*asin(1.)
91  avocado = 6.02214179e23   ! added by RW
[1672]92  kbol = 1.38064852e-23  ! added by JVO for Titan chem
[135]93
[1524]94  ! Initialize some "temporal and calendar" related variables
[2366]95#ifndef MESOSCALE
[1525]96  CALL init_time(day_ini,pdaysec,nday,ptimestep)
[2366]97#endif
[135]98
[1525]99  ! read in some parameters from "run.def" for physics,
100  ! or shared between dynamics and physics.
101  call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
102                                    ! in dynamical steps
103  call getin_p("day_step",day_step) ! number of dynamical steps per day
104  call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
[135]105
[1670]106  ! do we read a startphy.nc file? (default: .true.)
107  call getin_p("startphy_file",startphy_file)
108 
[135]109! --------------------------------------------------------------
110!  Reading the "callphys.def" file controlling some key options
111! --------------------------------------------------------------
112     
[1524]113  ! check that 'callphys.def' file is around
114  OPEN(99,file='callphys.def',status='old',form='formatted',iostat=ierr)
115  CLOSE(99)
116  IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module
[135]117     
[1315]118!!!      IF(ierr.EQ.0) THEN
[1524]119  IF(iscallphys) THEN
120     PRINT*
121     PRINT*
122     PRINT*,'--------------------------------------------'
123     PRINT*,' inifis: Parametres pour la physique (callphys.def)'
124     PRINT*,'--------------------------------------------'
[135]125
[1524]126     write(*,*) "Directory where external input files are:"
127     ! default 'datadir' is set in "datadir_mod"
128     call getin_p("datadir",datadir) ! default path
129     write(*,*) " datadir = ",trim(datadir)
[374]130
[1524]131     write(*,*) "Run with or without tracer transport ?"
132     tracer=.false. ! default value
133     call getin_p("tracer",tracer)
134     write(*,*) " tracer = ",tracer
[135]135
[1524]136     write(*,*) "Run with or without atm mass update ", &
137            " due to tracer evaporation/condensation?"
138     mass_redistrib=.false. ! default value
139     call getin_p("mass_redistrib",mass_redistrib)
140     write(*,*) " mass_redistrib = ",mass_redistrib
[728]141
[1524]142     write(*,*) "Diurnal cycle ?"
143     write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
144     diurnal=.true. ! default value
145     call getin_p("diurnal",diurnal)
146     write(*,*) " diurnal = ",diurnal
[135]147
[1524]148     write(*,*) "Seasonal cycle ?"
149     write(*,*) "(if season=false, Ls stays constant, to value ", &
150         "set in 'start'"
151     season=.true. ! default value
152     call getin_p("season",season)
153     write(*,*) " season = ",season
[2366]154     
155     write(*,*) "No seasonal cycle: initial day to lock the run during restart"
156     noseason_day=0.0 ! default value
157     call getin_p("noseason_day",noseason_day)
158     write(*,*) "noseason_day=", noseason_day
[135]159
[1524]160     write(*,*) "Tidally resonant rotation ?"
161     tlocked=.false. ! default value
162     call getin_p("tlocked",tlocked)
163     write(*,*) "tlocked = ",tlocked
[135]164
[1524]165     write(*,*) "Saturn ring shadowing ?"
166     rings_shadow = .false.
167     call getin_p("rings_shadow", rings_shadow)
168     write(*,*) "rings_shadow = ", rings_shadow
[1133]169         
[1524]170     write(*,*) "Compute latitude-dependent gravity field?"
171     oblate = .false.
172     call getin_p("oblate", oblate)
173     write(*,*) "oblate = ", oblate
[1194]174
[1524]175     write(*,*) "Flattening of the planet (a-b)/a "
176     flatten = 0.0
177     call getin_p("flatten", flatten)
[1672]178     write(*,*) "flatten = ", flatten         
[1194]179
[1524]180     write(*,*) "Needed if oblate=.true.: J2"
181     J2 = 0.0
182     call getin_p("J2", J2)
183     write(*,*) "J2 = ", J2
[1194]184         
[1524]185     write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
186     MassPlanet = 0.0
187     call getin_p("MassPlanet", MassPlanet)
188     write(*,*) "MassPlanet = ", MassPlanet         
[1194]189
[1524]190     write(*,*) "Needed if oblate=.true.: Planet mean radius (m)"
191     Rmean = 0.0
192     call getin_p("Rmean", Rmean)
193     write(*,*) "Rmean = ", Rmean
[1947]194     
195     write(*,*) "Compute effective altitude-dependent gravity field?"
196     eff_gz = .false.
197     call getin_p("eff_gz", eff_gz)
198     write(*,*) "eff_gz = ", eff_gz
[3083]199
200     write(*,*) "Nudging of zonal wind to keep the super-rotation"
201     nudging_u = .true.
202     call getin_p("nudging_u", nudging_u)
203     write(*,*) "nudging_u = ", nudging_u
204
205     write(*,*) "Nudging relaxation time (here ~ 1 titanian year in secondes)"
206     nudging_dt = 9.3e8
207     call getin_p("nudging_dt", nudging_dt)
208     write(*,*) "nudging_dt = ", nudging_dt
209
210     write(*,*) "New optics for clouds"
211     opt4clouds = .true.
212     call getin_p("opt4clouds", opt4clouds)
213     write(*,*) "opt4clouds = ", opt4clouds
214
215     write(*,*) "Cloudy fraction in the Clear / Dark column method"
216     Fcloudy = 1.0
217     call getin_p("Fcloudy", Fcloudy)
218     write(*,*) "Fcloudy = ", Fcloudy
219
[3318]220     write(*,*) "Adjustment of droplet properties in the VIS"
221     FCVIS = 1.0
222     call getin_p("FCVIS", FCVIS)
223     write(*,*) "FCVIS = ", FCVIS
224
225     write(*,*) "Adjustment of droplet properties in the IR"
226     FCIR = 1.0
227     call getin_p("FCIR", FCIR)
228     write(*,*) "FCIR = ", FCIR
[2241]229     
[3083]230     write(*,*) "Adjustment of aerosol properties in the VIS"
231     FHVIS = 1.0
232     call getin_p("FHVIS", FHVIS)
233     write(*,*) "FHVIS = ", FHVIS
234
235     write(*,*) "Adjustment of aerosol properties in the IR"
236     FHIR = 1.0
237     call getin_p("FHIR", FHIR)
238     write(*,*) "FHIR = ", FHIR
239
240     write(*,*) "Infinite tank of CH4 ?"
241     resCH4_inf = .false.
242     call getin_p("resCH4_inf", resCH4_inf)
243     write(*,*) "resCH4_inf = ", resCH4_inf
244     
[2241]245     ! sanity check warning
246     if (eff_gz) then
247       print*,"WARNING : You run chemistry with effective altitude-dependent gravity field !!"
248       print*,"You will have no coherence in your heating rates between physics and dynamics !!"
249       print*,"I let you continue but you should rather set eff_gz =.false. ..."
250     endif
251
[1194]252         
[135]253! Test of incompatibility:
254! if tlocked, then diurnal should be false
[1524]255     if (tlocked.and.diurnal) then
256       print*,'If diurnal=true, we should turn off tlocked.'
257       stop
258     endif
[135]259
[1524]260     write(*,*) "Tidal resonance ratio ?"
261     nres=0          ! default value
262     call getin_p("nres",nres)
263     write(*,*) "nres = ",nres
[135]264
[1524]265     write(*,*) "Write some extra output to the screen ?"
266     lwrite=.false. ! default value
267     call getin_p("lwrite",lwrite)
268     write(*,*) " lwrite = ",lwrite
[135]269
[1524]270     write(*,*) "Save statistics in file stats.nc ?"
271     callstats=.true. ! default value
272     call getin_p("callstats",callstats)
273     write(*,*) " callstats = ",callstats
[135]274
[1524]275     write(*,*) "Test energy conservation of model physics ?"
276     enertest=.false. ! default value
277     call getin_p("enertest",enertest)
278     write(*,*) " enertest = ",enertest
[135]279
[1524]280     write(*,*) "Check to see if cpp values used match gases.def ?"
281     check_cpp_match=.true. ! default value
282     call getin_p("check_cpp_match",check_cpp_match)
283     write(*,*) " check_cpp_match = ",check_cpp_match
[538]284
[1524]285     write(*,*) "call radiative transfer ?"
286     callrad=.true. ! default value
287     call getin_p("callrad",callrad)
288     write(*,*) " callrad = ",callrad
[135]289
[1524]290     write(*,*) "call correlated-k radiative transfer ?"
291     corrk=.true. ! default value
292     call getin_p("corrk",corrk)
293     write(*,*) " corrk = ",corrk
[1822]294     
295     if (corrk) then
296       ! default path is set in datadir
297       write(*,*) "callcorrk: Correlated-k data base folder:",trim(datadir)
298       call getin_p("corrkdir",corrkdir)
299       write(*,*) " corrkdir = ",corrkdir
[2050]300       
301       write(*,*) "use correlated-k recombination instead of pre-mixed values ?"
[2241]302       corrk_recombin=.false.! default value
[2050]303       call getin_p("corrk_recombin",corrk_recombin)
304       write(*,*) " corrk_recombin = ",corrk_recombin
[1822]305     endif
306     
307     if (corrk .and. ngrid.eq.1) then
308       write(*,*) "simulate global averaged conditions ?"
309       global1d = .false. ! default value
310       call getin_p("global1d",global1d)
311       write(*,*) " global1d = ",global1d
312       
313       ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
314       if (global1d.and.diurnal) then
315          write(*,*) "if global1d is true, diurnal must be set to false"
316          stop
317       endif
[135]318
[1822]319       if (global1d) then
320          write(*,*) "Solar Zenith angle (deg.) ?"
321          write(*,*) "(assumed for averaged solar flux S/4)"
322          szangle=60.0  ! default value
323          call getin_p("szangle",szangle)
324          write(*,*) " szangle = ",szangle
325       endif
326     endif
327
[1524]328     write(*,*) "prohibit calculations outside corrk T grid?"
329     strictboundcorrk=.true. ! default value
330     call getin_p("strictboundcorrk",strictboundcorrk)
331     write(*,*) "strictboundcorrk = ",strictboundcorrk
[1145]332
[2366]333     write(*,*) "Minimum atmospheric temperature for Planck function integration ?"
334     tplanckmin=30.0 ! default value
335     call getin_p("tplanckmin",tplanckmin)
336     write(*,*) " tplanckmin = ",tplanckmin
337 
338     write(*,*) "Maximum atmospheric temperature for Planck function integration ?"
339     tplanckmax=1500.0 ! default value
340     call getin_p("tplanckmax",tplanckmax)
341     write(*,*) " tplanckmax = ",tplanckmax
342 
343     write(*,*) "Temperature step for Planck function integration ?"
344     dtplanck=0.1 ! default value
345     call getin_p("dtplanck",dtplanck)
346     write(*,*) " dtplanck = ",dtplanck
347 
[1524]348     write(*,*) "call gaseous absorption in the visible bands?", &
349                    "(matters only if callrad=T)"
350     callgasvis=.false. ! default value
351     call getin_p("callgasvis",callgasvis)
352     write(*,*) " callgasvis = ",callgasvis
[538]353       
[1524]354     write(*,*) "call continuum opacities in radiative transfer ?", &
355                    "(matters only if callrad=T)"
356     continuum=.true. ! default value
357     call getin_p("continuum",continuum)
358     write(*,*) " continuum = ",continuum
[538]359 
[2245]360     write(*,*) "version for H2H2 CIA file ?"
361     versH2H2cia=2011 ! default value (should be 2018 but retrocompatibility first)
362     call getin_p("versH2H2cia",versH2H2cia)
363     write(*,*) " versH2H2cia = ",versH2H2cia
364     ! Sanity check
365     if (versH2H2cia.ne.2011 .and. versH2H2cia.ne.2018) then
366        print*,'Error: Choose a correct value (2011 or 2018) for versH2H2cia !'
367        call abort
368     endif
369
[1524]370     write(*,*) "call turbulent vertical diffusion ?"
371     calldifv=.true. ! default value
372     call getin_p("calldifv",calldifv)
373     write(*,*) " calldifv = ",calldifv
[135]374
[1524]375     write(*,*) "use turbdiff instead of vdifc ?"
376     UseTurbDiff=.true. ! default value
377     call getin_p("UseTurbDiff",UseTurbDiff)
378     write(*,*) " UseTurbDiff = ",UseTurbDiff
[596]379
[1524]380     write(*,*) "call convective adjustment ?"
381     calladj=.true. ! default value
382     call getin_p("calladj",calladj)
383     write(*,*) " calladj = ",calladj
[135]384
[1524]385     write(*,*) "Radiative timescale for Newtonian cooling ?"
386     tau_relax=30. ! default value
387     call getin_p("tau_relax",tau_relax)
388     write(*,*) " tau_relax = ",tau_relax
389     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
[253]390
[1524]391     write(*,*)"call thermal conduction in the soil ?"
392     callsoil=.true. ! default value
393     call getin_p("callsoil",callsoil)
394     write(*,*) " callsoil = ",callsoil
[135]395         
[1524]396     write(*,*)"Rad transfer is computed every iradia", &
397                   " physical timestep"
398     iradia=1 ! default value
399     call getin_p("iradia",iradia)
400     write(*,*)" iradia = ",iradia
[135]401       
[1524]402     write(*,*)"Rayleigh scattering ?"
403     rayleigh=.false.
404     call getin_p("rayleigh",rayleigh)
405     write(*,*)" rayleigh = ",rayleigh
[135]406
[1524]407     write(*,*) "Use blackbody for stellar spectrum ?"
408     stelbbody=.false. ! default value
409     call getin_p("stelbbody",stelbbody)
410     write(*,*) " stelbbody = ",stelbbody
[135]411
[1524]412     write(*,*) "Stellar blackbody temperature ?"
413     stelTbb=5800.0 ! default value
414     call getin_p("stelTbb",stelTbb)
415     write(*,*) " stelTbb = ",stelTbb
[253]416
[1524]417     write(*,*)"Output mean OLR in 1D?"
418     meanOLR=.false.
419     call getin_p("meanOLR",meanOLR)
420     write(*,*)" meanOLR = ",meanOLR
[135]421
[1524]422     write(*,*)"Output spectral OLR in 3D?"
423     specOLR=.false.
424     call getin_p("specOLR",specOLR)
425     write(*,*)" specOLR = ",specOLR
[135]426
[2138]427     write(*,*)"Output diagnostic optical thickness attenuation exp(-tau)?"
[2131]428     diagdtau=.false.
429     call getin_p("diagdtau",diagdtau)
430     write(*,*)" diagdtau = ",diagdtau
431
[1524]432     write(*,*)"Uniform absorption in radiative transfer?"
433     graybody=.false.
434     call getin_p("graybody",graybody)
435     write(*,*)" graybody = ",graybody
[253]436
[1672]437! Chemistry
438
439     write(*,*) "Run with or without chemistry?"
440     callchim=.false. ! default value
441     call getin_p("callchim",callchim)
442     write(*,*) " callchim = ",callchim
443
444     ! sanity check
445     if (callchim.and.(.not.tracer)) then
446       print*,"You are running chemistry without tracer"
447       print*,"Please start again with tracer =.true."
448       stop
449     endif
[1947]450     
[1672]451     write(*,*)"Chemistry is computed every ichim", &
452                   " physical timestep"
453     ichim=1 ! default value
454     call getin_p("ichim",ichim)
455     write(*,*)" ichim = ",ichim
456
[1795]457! Microphysics
458
[2046]459     write(*,*) "Use haze seasonal model from Karkoschka 2016 ?"
460     seashaze=.true. ! default value
461     call getin_p("seashaze",seashaze)
462     write(*,*)" seashaze = ",seashaze
463
[1795]464     write(*,*) "Run with or without microphysics?"
465     callmufi=.false. ! default value
466     call getin_p("callmufi",callmufi)
467     write(*,*)" callmufi = ",callmufi
468
469     ! sanity check
470     if (callmufi.and.(.not.tracer)) then
471       print*,"You are running microphysics without tracer"
472       print*,"Please start again with tracer =.true."
473       stop
474     endif
475
476     write(*,*) "Compute clouds?"
477     callclouds=.false. ! default value
478     call getin_p("callclouds",callclouds)
479     write(*,*)" callclouds = ",callclouds
480
481     ! sanity check
482     if (callclouds.and.(.not.callmufi)) then
483       print*,"You are trying to make clouds without microphysics !"
484       print*,"Please start again with callmufi =.true."
485       stop
486     endif
487
[1897]488     write(*,*) "Disable the coupling of microphysics within rad. transf. ?"
489     write(*,*) "If disabled we will assume a planetwide vert. profile of extinction ..."
490     uncoupl_optic_haze=.true. ! default value - true as long as the microphysics is bugged
491     call getin_p("uncoupl_optic_haze",uncoupl_optic_haze)
492     write(*,*)" uncoupl_optic_haze = ",uncoupl_optic_haze
[1795]493
494     write(*,*) "Pressure level of aer. production (Pa) ?"
495     p_prod=1.0 ! default value
496     call getin_p("p_prod",p_prod)
[1822]497     write(*,*)" p_prod = ",p_prod
[1795]498     
499     write(*,*) "Aerosol production rate (kg.m-2.s-1) ?"
500     tx_prod=3.5e-13 ! default value
501     call getin_p("tx_prod",tx_prod)
[1822]502     write(*,*)" tx_prod = ",tx_prod
[1795]503
504     write(*,*) "Equivalent radius production (m) ?"
505     rc_prod=2.0e-8 ! default value
506     call getin_p("rc_prod",rc_prod)
[1822]507     write(*,*)" rhc_prod = ",rc_prod
[1795]508
509     write(*,*) "Radius of air (nitrogen) molecule (m) ?"
510     air_rad=1.75e-10 ! default value
511     call getin_p("air_rad",air_rad)
[1822]512     write(*,*)" air_rad = ",air_rad
[1795]513
514     write(*,*) "Path to microphys. config file ?"
515     config_mufi='datagcm/microphysics/config.cfg' ! default value
516     call getin_p("config_mufi",config_mufi)
[1822]517     write(*,*)" config_mufi = ",config_mufi
[1795]518
[1538]519! Soil model
520     write(*,*)"Number of sub-surface layers for soil scheme?"
521     ! default value of nsoilmx set in comsoil_h
522     call getin_p("nsoilmx",nsoilmx)
523     write(*,*)" nsoilmx=",nsoilmx
524     
525     write(*,*)"Thickness of topmost soil layer (m)?"
526     ! default value of lay1_soil set in comsoil_h
527     call getin_p("lay1_soil",lay1_soil)
528     write(*,*)" lay1_soil = ",lay1_soil
529     
530     write(*,*)"Coefficient for soil layer thickness distribution?"
531     ! default value of alpha_soil set in comsoil_h
532     call getin_p("alpha_soil",alpha_soil)
533     write(*,*)" alpha_soil = ",alpha_soil
534
[1524]535     write(*,*)"Remove lower boundary?"
536     nosurf=.false.
537     call getin_p("nosurf",nosurf)
538     write(*,*)" nosurf = ",nosurf
[253]539
540! Tests of incompatibility:
[1524]541     if (nosurf.and.callsoil) then
542       print*,'nosurf not compatible with soil scheme!'
543       print*,'... got to make a choice!'
544       call abort
545     endif
[253]546
[1524]547     write(*,*)"Add an internal heat flux?", &
548                   "... matters only if callsoil=F"
549     intheat=0.
550     call getin_p("intheat",intheat)
551     write(*,*)" intheat = ",intheat
[952]552
[1524]553     write(*,*)"Use Newtonian cooling for radiative transfer?"
554     newtonian=.false.
555     call getin_p("newtonian",newtonian)
556     write(*,*)" newtonian = ",newtonian
[253]557
558! Tests of incompatibility:
[1524]559     if (newtonian.and.corrk) then
560       print*,'newtonian not compatible with correlated-k!'
561       call abort
562     endif
563     if (newtonian.and.calladj) then
564       print*,'newtonian not compatible with adjustment!'
565       call abort
566     endif
567     if (newtonian.and.calldifv) then
568       print*,'newtonian not compatible with a boundary layer!'
569       call abort
570     endif
[253]571
[1524]572     write(*,*)"Test physics timescale in 1D?"
573     testradtimes=.false.
574     call getin_p("testradtimes",testradtimes)
575     write(*,*)" testradtimes = ",testradtimes
[253]576
577! Test of incompatibility:
578! if testradtimes used, we must be in 1D
[1524]579     if (testradtimes.and.(ngrid.gt.1)) then
580       print*,'testradtimes can only be used in 1D!'
581       call abort
582     endif
[253]583
[1524]584     write(*,*)"Which star?"
585     startype=1 ! default value = Sol
586     call getin_p("startype",startype)
587     write(*,*)" startype = ",startype
[135]588
[1524]589     write(*,*)"Value of stellar flux at 1 AU?"
590     Fat1AU=1356.0 ! default value = Sol today
591     call getin_p("Fat1AU",Fat1AU)
592     write(*,*)" Fat1AU = ",Fat1AU
[135]593
[1524]594     write(*,*) "Does user want to force cpp and mugaz?"
595     force_cpp=.false. ! default value
596     call getin_p("force_cpp",force_cpp)
597     write(*,*) " force_cpp = ",force_cpp
[589]598
[1524]599     if (force_cpp) then
600       mugaz = -99999.
601       PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
602       call getin_p("mugaz",mugaz)
603       IF (mugaz.eq.-99999.) THEN
604           PRINT *, "mugaz must be set if force_cpp = T"
605           STOP
606       ELSE
607           write(*,*) "mugaz=",mugaz
608       ENDIF
609       cpp = -99999.
610       PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
611       call getin_p("cpp",cpp)
612       IF (cpp.eq.-99999.) THEN
613           PRINT *, "cpp must be set if force_cpp = T"
614           STOP
615       ELSE
616           write(*,*) "cpp=",cpp
617       ENDIF
618     endif ! of if (force_cpp)
[1648]619     
620     call su_gases(nlayer,tracer)     
[2366]621     call calc_cpp_mugaz
622
[1524]623     PRINT*,'--------------------------------------------'
624     PRINT*
625     PRINT*
626  ELSE
627     write(*,*)
628     write(*,*) 'Cannot read file callphys.def. Is it here ?'
629     stop
630  ENDIF ! of IF(iscallphys)
[135]631
[1524]632  PRINT*
633  PRINT*,'inifis: daysec',daysec
634  PRINT*
635  PRINT*,'inifis: The radiative transfer is computed:'
636  PRINT*,'           each ',iradia,' physical time-step'
637  PRINT*,'        or each ',iradia*dtphys,' seconds'
638  PRINT*
[1672]639  PRINT*,'inifis: Chemistry is computed:'
640  PRINT*,'           each ',ichim,' physical time-step'
641  PRINT*,'        or each ',ichim*dtphys,' seconds'
642  PRINT*
[135]643
644!-----------------------------------------------------------------------
645!     Some more initialization:
646!     ------------------------
647
[1542]648  ! Initializations for comgeomfi_h
[2366]649#ifndef MESOSCALE
[1542]650  totarea=SSUM(ngrid,parea,1)
651  call planetwide_sumval(parea,totarea_planet)
[787]652
[1524]653  !! those are defined in comdiurn_h.F90
654  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
655  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
656  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
657  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
[787]658
[1524]659  DO ig=1,ngrid
660     sinlat(ig)=sin(plat(ig))
661     coslat(ig)=cos(plat(ig))
662     sinlon(ig)=sin(plon(ig))
663     coslon(ig)=cos(plon(ig))
664  ENDDO
[2366]665#endif
[1722]666  ! initialize variables in radinc_h
[2366]667  call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
[1529]668 
[1524]669  ! allocate "comsoil_h" arrays
670  call ini_comsoil_h(ngrid)
[2366]671   
[1524]672  END SUBROUTINE inifis
[135]673
[1524]674END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.