source: trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 @ 3914

Last change on this file since 3914 was 3894, checked in by debatzbr, 4 months ago

Pluto PCM: Allows for haze production with vmr CH4 fix profile.
BBT

File size: 54.4 KB
Line 
1MODULE inifis_mod
2IMPLICIT NONE
3
4CONTAINS
5
6  SUBROUTINE inifis(ngrid,nlayer,nq, &
7             day_ini,pdaysec,nday,ptimestep, &
8             plat,plon,parea, &
9             prad,pg,pr,pcpp)
10
11  use init_print_control_mod, only: init_print_control
12  use radinc_h, only: ini_radinc_h, naerkind
13  use radcommon_h, only: ini_radcommon_h
14  use radii_mod, only: radfixed
15  use datafile_mod, only: datadir,hazeprop_file,hazerad_file,hazemmr_file,hazedens_file, &
16                          config_mufi, aersprop_file, aerfprop_file
17  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
18  use comgeomfi_h, only: totarea, totarea_planet
19  use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
20  use time_phylmdz_mod, only: diagfi_output_rate,startfi_output_rate, &
21                              init_time, daysec, dtphys
22  use comcstfi_mod, only: rad, cpp, g, r, rcp, &
23                          mugaz, pi, avocado
24  use planetwide_mod, only: planetwide_sumval
25  use callkeys_mod
26  use surfdat_h
27  use wstats_mod, only: callstats
28  use writediagsoil_mod, only: diagsoil
29  use ioipsl_getin_p_mod, only : getin_p
30  use mod_phys_lmdz_para, only : is_parallel, is_master, bcast
31
32!=======================================================================
33!
34!   purpose:
35!   -------
36!
37!   Initialisation for the physical parametrisations of the LMD
38!   Generic Model.
39!
40!   author: Frederic Hourdin 15 / 10 /93
41!   -------
42!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
43!             Ehouarn Millour (oct. 2008) tracers are now identified
44!              by their names and may not be contiguously
45!              stored in the q(:,:,:,:) array
46!             E.M. (june 2009) use getin routine to load parameters
47!
48!
49!   arguments:
50!   ----------
51!
52!   input:
53!   ------
54!
55!    ngrid                 Size of the horizontal grid.
56!                          All internal loops are performed on that grid.
57!    nlayer                Number of vertical layers.
58!    pdayref               Day of reference for the simulation
59!    pday                  Number of days counted from the North. Spring
60!                          equinoxe.
61!
62!=======================================================================
63!
64!-----------------------------------------------------------------------
65!   declarations:
66!   -------------
67  use datafile_mod, only: datadir
68  use ioipsl_getin_p_mod, only: getin_p
69  IMPLICIT NONE
70
71
72
73  REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
74  INTEGER,INTENT(IN) :: nday
75  INTEGER,INTENT(IN) :: ngrid,nlayer,nq
76  REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
77  integer,intent(in) :: day_ini
78
79  INTEGER ig
80  CHARACTER(len=20) :: rname="inifis" ! routine name, for messages
81
82  EXTERNAL iniorbit,orbite
83  EXTERNAL SSUM
84  REAL SSUM
85
86  ! deprecated parameters
87  REAL :: ecritphy ! to check that this obsolete flag is no longer used...
88  logical aerohaze
89
90  ! Initialize flags lunout, prt_level, debug (in print_control_mod)
91  CALL init_print_control
92
93  ! initialize constants in comcstfi_mod
94  rad=prad
95  cpp=pcpp
96  g=pg
97  r=pr
98  rcp=r/cpp
99  mugaz=8.314*1000./pr ! dummy init
100  pi=2.*asin(1.)
101  avocado = 6.02214179e23   ! added by RW
102
103  ! Initialize some "temporal and calendar" related variables
104#ifndef MESOSCALE
105  CALL init_time(day_ini,pdaysec,nday,ptimestep)
106#endif
107
108
109  ! read in some parameters from "run.def" for physics,
110  ! or shared between dynamics and physics.
111  ecritphy=-666 ! dummy default value
112  call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
113                                    ! in dynamical steps
114  if (ecritphy/=-666) then
115    call abort_physic(rname, &
116         "Error: parameter ecritphy is obsolete! Remove it! "//&
117         "And use diagfi_output_rate instead",1)
118  endif
119
120  ! do we read a startphy.nc file? (default: .true.)
121  call getin_p("startphy_file",startphy_file)
122
123! --------------------------------------------------------------
124!  Reading the "callphys.def" file controlling some key options
125! --------------------------------------------------------------
126
127  IF (is_master) THEN
128    ! check if 'callphys.def' file is around
129    inquire(file="callphys.def",exist=iscallphys)
130    write(*,*) trim(rname)//": iscallphys=",iscallphys
131  ENDIF
132  call bcast(iscallphys)
133
134  IF(iscallphys) THEN
135
136     if (is_master) then
137       write(*,*)
138       write(*,*)'--------------------------------------------'
139       write(*,*)trim(rname)//': Parameters for the physics (callphys.def)'
140       write(*,*)'--------------------------------------------'
141     endif
142
143     if (is_master) write(*,*) trim(rname)//&
144       ": Directory where external input files are:"
145     ! default 'datadir' is set in "datadir_mod"
146     call getin_p("datadir",datadir) ! default path
147     if (is_master) write(*,*) trim(rname)//": datadir = ",trim(datadir)
148
149     if (is_master) write(*,*) "Haze optical properties datafile"
150     hazeprop_file="optprop_tholins_fractal050nm.dat"  ! default file
151     call getin_p("hazeprop_file",hazeprop_file)
152     if (is_master) write(*,*) trim(rname)//" hazeprop_file = ",trim(hazeprop_file)
153
154     if (is_master) write(*,*) "Haze radii datafile"
155     hazerad_file="hazerad.txt"  ! default file
156     call getin_p("hazerad_file",hazerad_file)
157     if (is_master) write(*,*) trim(rname)//" hazerad_file = ",trim(hazerad_file)
158
159     if (is_master) write(*,*) "Haze mmr datafile"
160     hazemmr_file="None"  ! default file
161     call getin_p("hazemmr_file",hazemmr_file)
162     if (is_master) write(*,*) trim(rname)//" hazemmr_file = ",trim(hazemmr_file)
163     if (is_master) write(*,*) "Haze density datafile"
164     hazedens_file="None"  ! default file
165     call getin_p("hazedens_file",hazedens_file)
166     if (is_master) write(*,*) trim(rname)//" hazedens_file = ",trim(hazedens_file)
167
168     if (is_master) write(*,*) trim(rname)//&
169       ": Run with or without tracer transport ?"
170     tracer=.false. ! default value
171     call getin_p("tracer",tracer)
172     if (is_master) write(*,*) trim(rname)//": tracer = ",tracer
173
174     if (is_master) write(*,*) trim(rname)//&
175       ": Output rate for diagfi.nc file (in physics steps) ?"
176     diagfi_output_rate=24 !default value
177     call getin_p("diagfi_output_rate",diagfi_output_rate)
178     if (is_master) write(*,*) trim(rname)//": diagfi_output_rate = ",&
179                               diagfi_output_rate
180
181     if (is_master) write(*,*) trim(rname)//&
182       ": Output rate for start/startfi.nc files (in physics steps) ?"
183     startfi_output_rate=0 !default value
184     call getin_p("startfi_output_rate",startfi_output_rate)
185     if (is_master) write(*,*) trim(rname)//": startfi_output_rate = ",&
186                               startfi_output_rate
187
188     if (is_master) write(*,*) trim(rname)//&
189       ": Run with or without atm mass update "//&
190       " due to tracer evaporation/condensation?"
191     mass_redistrib=.false. ! default value
192     call getin_p("mass_redistrib",mass_redistrib)
193     if (is_master) write(*,*) trim(rname)//&
194       ": mass_redistrib = ",mass_redistrib
195
196     if (is_master) then
197       write(*,*) trim(rname)//": Diurnal cycle ?"
198       write(*,*) trim(rname)//&
199       ": (if diurnal=false, diurnal averaged solar heating)"
200     endif
201     diurnal=.true. ! default value
202     call getin_p("diurnal",diurnal)
203     if (is_master) write(*,*) trim(rname)//": diurnal = ",diurnal
204
205     if (is_master) then
206       write(*,*) trim(rname)//": Seasonal cycle ?"
207       write(*,*) trim(rname)//&
208       ": (if season=false, Ls stays constant, to value ", &
209       "set in 'start'"
210     endif
211     season=.true. ! default value
212     call getin_p("season",season)
213     if (is_master) write(*,*) trim(rname)//": season = ",season
214
215     if (is_master) then
216       write(*,*) trim(rname)//": Fast mode (nogcm) ?"
217     endif
218     fast=.false. ! default value
219     call getin_p("fast",fast)
220     if (is_master) write(*,*) trim(rname)//": fast = ",fast
221
222     if (is_master) write(*,*) trim(rname)//"Run with Triton orbit ?"
223     triton=.false. ! default value
224     call getin_p("triton",triton)
225     if (is_master) write(*,*) trim(rname)//" triton = ",triton
226     convergeps=.false. ! default value
227     call getin_p("convergeps",convergeps)
228     if (is_master) write(*,*) trim(rname)//" convergeps = ",convergeps
229
230     conservn2=.false. ! default value
231     call getin_p("conservn2",conservn2)
232     if (is_master) write(*,*) trim(rname)//" conservn2 = ",conservn2
233     conservch4=.false. ! default value
234     call getin_p("conservch4",conservch4)
235     if (is_master) write(*,*) trim(rname)//" conservch4 = ",conservch4
236
237     if (is_master) write(*,*) trim(rname)//"KBO runs (eris, makemake) ?"
238     kbo=.false. ! default value
239     call getin_p("kbo",kbo)
240     if (is_master) write(*,*) trim(rname)//" kbo = ",kbo
241
242     if (is_master) write(*,*) trim(rname)//"Specific paleo run ?"
243     paleo=.false. ! default value
244     call getin_p("paleo",paleo)
245     if (is_master) write(*,*) trim(rname)//" paleo = ",paleo
246
247     if (is_master) write(*,*) trim(rname)//"paleoclimate step (Earth years) "
248     paleoyears=10000. ! default value
249     call getin_p("paleoyears",paleoyears)
250     if (is_master) write(*,*) trim(rname)//" paleoyears = ",paleoyears
251
252     if (is_master) write(*,*) trim(rname)//&
253       ": No seasonal cycle: initial day to lock the run during restart"
254     noseason_day=0.0 ! default value
255     call getin_p("noseason_day",noseason_day)
256     if (is_master) write(*,*) trim(rname)//": noseason_day=", noseason_day
257
258     if (is_master) write(*,*) trim(rname)//&
259       ": Write some extra output to the screen ?"
260     lwrite=.false. ! default value
261     call getin_p("lwrite",lwrite)
262     if (is_master) write(*,*) trim(rname)//": lwrite = ",lwrite
263
264     if (is_master) write(*,*) trim(rname)//&
265       ": Save statistics in file stats.nc ?"
266     callstats=.true. ! default value
267     call getin_p("callstats",callstats)
268     if (is_master) write(*,*) trim(rname)//": callstats = ",callstats
269
270     if (is_master) write(*,*) trim(rname)//&
271       ": Write sub-surface fields in file diagsoil.nc ?"
272     diagsoil=.false. ! default value
273     call getin_p("diagsoil",diagsoil)
274     if (is_master) write(*,*) trim(rname)//" diagsoil = ",diagsoil
275
276     if (is_master) write(*,*) trim(rname)//&
277       ": Test energy conservation of model physics ?"
278     enertest=.false. ! default value
279     call getin_p("enertest",enertest)
280     if (is_master) write(*,*) trim(rname)//": enertest = ",enertest
281
282     if (is_master) write(*,*) trim(rname)//": call radiative transfer ?"
283     callrad=.true. ! default value
284     call getin_p("callrad",callrad)
285     if (is_master) write(*,*) trim(rname)//": callrad = ",callrad
286
287     if (is_master) write(*,*) trim(rname)//&
288       ": call correlated-k radiative transfer ?"
289     corrk=.true. ! default value
290     call getin_p("corrk",corrk)
291     if (is_master) write(*,*) trim(rname)//": corrk = ",corrk
292
293     if (is_master) write(*,*) trim(rname)//&
294       ": prohibit calculations outside corrk T grid?"
295     strictboundcorrk=.true. ! default value
296     call getin_p("strictboundcorrk",strictboundcorrk)
297     if (is_master) write(*,*) trim(rname)//&
298       ": strictboundcorrk = ",strictboundcorrk
299
300     if (is_master) write(*,*) trim(rname)//&
301       ": prohibit calculations outside CIA T grid?"
302     strictboundcia=.true. ! default value
303     call getin_p("strictboundcia",strictboundcia)
304     if (is_master) write(*,*) trim(rname)//&
305       ": strictboundcia = ",strictboundcia
306
307     if (is_master) write(*,*) trim(rname)//&
308       ": Minimum atmospheric temperature for Planck function integration ?"
309     tplanckmin=30.0 ! default value
310     call getin_p("tplanckmin",tplanckmin)
311     if (is_master) write(*,*) trim(rname)//": tplanckmin = ",tplanckmin
312
313     if (is_master) write(*,*) trim(rname)//&
314       ": Maximum atmospheric temperature for Planck function integration ?"
315     tplanckmax=1500.0 ! default value
316     call getin_p("tplanckmax",tplanckmax)
317     if (is_master) write(*,*) trim(rname)//": tplanckmax = ",tplanckmax
318
319     if (is_master) write(*,*) trim(rname)//&
320       ": Temperature step for Planck function integration ?"
321     dtplanck=0.1 ! default value
322     call getin_p("dtplanck",dtplanck)
323     if (is_master) write(*,*) trim(rname)//": dtplanck = ",dtplanck
324
325     if (is_master) write(*,*) trim(rname)//&
326       ": call gaseous absorption in the visible bands?"//&
327       " (matters only if callrad=T)"
328     callgasvis=.false. ! default value
329     call getin_p("callgasvis",callgasvis)
330     if (is_master) write(*,*) trim(rname)//": callgasvis = ",callgasvis
331
332     if (is_master) write(*,*) trim(rname)//&
333       ": call continuum opacities in radiative transfer ?"//&
334       " (matters only if callrad=T)"
335     continuum=.true. ! default value
336     call getin_p("continuum",continuum)
337     if (is_master) write(*,*) trim(rname)//": continuum = ",continuum
338
339     if (is_master) write(*,*) trim(rname)//&
340       ": call turbulent vertical diffusion ?"
341     calldifv=.false. ! default value
342     call getin_p("calldifv",calldifv)
343     if (is_master) write(*,*) trim(rname)//": calldifv = ",calldifv
344     vertdiff=.true. ! default value
345     call getin_p("vertdiff",vertdiff)
346     if (is_master) write(*,*) trim(rname)//": vertdiff = ",vertdiff
347
348     if (is_master) write(*,*) trim(rname)//": use turbdiff instead of vdifc ?"
349     UseTurbDiff=.false. ! default value
350     call getin_p("UseTurbDiff",UseTurbDiff)
351     if (is_master) write(*,*) trim(rname)//": UseTurbDiff = ",UseTurbDiff
352
353     if (is_master) write(*,*) trim(rname)//": call convective adjustment ?"
354     calladj=.false. ! default value
355     call getin_p("calladj",calladj)
356     if (is_master) write(*,*) trim(rname)//": calladj = ",calladj
357
358     if (is_master) write(*,*) trim(rname)//&
359     ": Radiative timescale for Newtonian cooling (in Earth days)?"
360     tau_relax=30. ! default value
361     call getin_p("tau_relax",tau_relax)
362     if (is_master) write(*,*) trim(rname)//": tau_relax = ",tau_relax
363     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
364
365     if (is_master) write(*,*)trim(rname)//&
366       ": call thermal conduction in the soil ?"
367     callsoil=.true. ! default value
368     call getin_p("callsoil",callsoil)
369     if (is_master) write(*,*) trim(rname)//": callsoil = ",callsoil
370
371     calllott=.false. ! default value
372     call getin_p("calllott",calllott)
373     if (is_master) write(*,*)trim(rname)//": calllott = ",calllott
374
375     if (is_master) write(*,*)trim(rname)//&
376       ": Rad transfer is computed every iradia", &
377       " physical timestep"
378     iradia=1 ! default value
379     call getin_p("iradia",iradia)
380     if (is_master) write(*,*)trim(rname)//": iradia = ",iradia
381
382     if (is_master) write(*,*)trim(rname)//": Rayleigh scattering ?"
383     rayleigh=.false.
384     call getin_p("rayleigh",rayleigh)
385     if (is_master) write(*,*)trim(rname)//": rayleigh = ",rayleigh
386
387     if (is_master) write(*,*) trim(rname)//&
388       ": Use blackbody for stellar spectrum ?"
389     stelbbody=.false. ! default value
390     call getin_p("stelbbody",stelbbody)
391     if (is_master) write(*,*) trim(rname)//": stelbbody = ",stelbbody
392
393     if (is_master) write(*,*) trim(rname)//": Stellar blackbody temperature ?"
394     stelTbb=5800.0 ! default value
395     call getin_p("stelTbb",stelTbb)
396     if (is_master) write(*,*) trim(rname)//": stelTbb = ",stelTbb
397
398     if (is_master) write(*,*) trim(rname)//": Output mean OLR in 1D?"
399     meanOLR=.false.
400     call getin_p("meanOLR",meanOLR)
401     if (is_master) write(*,*) trim(rname)//": meanOLR = ",meanOLR
402
403     if (is_master) write(*,*)trim(rname)//": Output spectral OLR in 3D?"
404     specOLR=.false.
405     call getin_p("specOLR",specOLR)
406     if (is_master) write(*,*)trim(rname)//": specOLR = ",specOLR
407
408     if (is_master) write(*,*)trim(rname)//": Operate in kastprof mode?"
409     kastprof=.false.
410     call getin_p("kastprof",kastprof)
411     if (is_master) write(*,*)trim(rname)//": kastprof = ",kastprof
412
413     if (is_master) write(*,*)trim(rname)//&
414       ": Uniform absorption in radiative transfer?"
415     graybody=.false.
416     call getin_p("graybody",graybody)
417     if (is_master) write(*,*)trim(rname)//": graybody = ",graybody
418
419! Soil model
420     if (is_master) write(*,*)trim(rname)//&
421       ": Number of sub-surface layers for soil scheme?"
422     ! default value of nsoilmx set in comsoil_h
423     call getin_p("nsoilmx",nsoilmx)
424     if (is_master) write(*,*)trim(rname)//": nsoilmx=",nsoilmx
425
426     if (is_master) write(*,*)trim(rname)//&
427       ": Thickness of topmost soil layer (m)?"
428     ! default value of lay1_soil set in comsoil_h
429     call getin_p("lay1_soil",lay1_soil)
430     if (is_master) write(*,*)trim(rname)//": lay1_soil = ",lay1_soil
431
432     if (is_master) write(*,*)trim(rname)//&
433       ": Coefficient for soil layer thickness distribution?"
434     ! default value of alpha_soil set in comsoil_h
435     call getin_p("alpha_soil",alpha_soil)
436     if (is_master) write(*,*)trim(rname)//": alpha_soil = ",alpha_soil
437
438! Chemistry in the thermosphere
439     if (is_master) write(*,*) trim(rname)//": Use deposition ?"
440     depos=.false.         ! default value
441     call getin_p("depos",depos)
442     if (is_master) write(*,*) trim(rname)//": depos = ",depos
443
444     if (is_master) write(*,*)trim(rname)//": Production of haze ?"
445     haze=.false. ! default value
446     call getin_p("haze",haze)
447     if (is_master) write(*,*)trim(rname)//": haze = ",haze
448
449
450
451      if (is_master) write(*,*)trim(rname)// "call thermal conduction ?"
452      callconduct=.false. ! default value
453      call getin_p("callconduct",callconduct)
454      if (is_master) write(*,*)trim(rname)// " callconduct = ",callconduct
455
456      if (is_master) write(*,*)trim(rname)// "call phitop ?"
457      phitop=0. ! default value
458      call getin_p("phitop",phitop)
459      if (is_master) write(*,*)trim(rname)// " phitop = ",phitop
460
461      if (is_master) write(*,*)trim(rname)// "call molecular viscosity ?"
462      callmolvis=.false. ! default value
463      call getin_p("callmolvis",callmolvis)
464      if (is_master) write(*,*)trim(rname)// " callmolvis = ",callmolvis
465
466      if (is_master) write(*,*)trim(rname)// "call molecular diffusion ?"
467      callmoldiff=.false. ! default value
468      call getin_p("callmoldiff",callmoldiff)
469      if (is_master) write(*,*)trim(rname)// " callmoldiff = ",callmoldiff
470
471! Global1D mean and solar zenith angle
472
473     if (ngrid.eq.1) then
474      PRINT*, 'Simulate global averaged conditions ?'
475      global1d = .false. ! default value
476      call getin_p("global1d",global1d)
477      write(*,*) "global1d = ",global1d
478
479      ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
480      if (global1d.and.diurnal) then
481         call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
482      endif
483
484      if (global1d) then
485         PRINT *,'Solar Zenith angle (deg.) ?'
486         PRINT *,'(assumed for averaged solar flux S/4)'
487         szangle=60.0  ! default value
488         call getin_p("szangle",szangle)
489         write(*,*) "szangle = ",szangle
490      endif
491   endif ! of if (ngrid.eq.1)
492
493! Test of incompatibility:
494! if kastprof used, we must be in 1D
495     if (kastprof.and.(ngrid.gt.1)) then
496       call abort_physic(rname,'kastprof can only be used in 1D!',1)
497     endif
498
499     if (is_master) write(*,*)trim(rname)//&
500       ": Stratospheric temperature for kastprof mode?"
501     Tstrat=167.0
502     call getin_p("Tstrat",Tstrat)
503     if (is_master) write(*,*)trim(rname)//": Tstrat = ",Tstrat
504
505     if (is_master) write(*,*)trim(rname)//": Remove lower boundary?"
506     nosurf=.false.
507     call getin_p("nosurf",nosurf)
508     if (is_master) write(*,*)trim(rname)//": nosurf = ",nosurf
509
510! Tests of incompatibility:
511     if (nosurf.and.callsoil) then
512       if (is_master) then
513         write(*,*)trim(rname)//'nosurf not compatible with soil scheme!'
514         write(*,*)trim(rname)//'... got to make a choice!'
515       endif
516       call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
517     endif
518
519     if (is_master) write(*,*)trim(rname)//": Add an internal heat flux?", &
520                   "... matters only if callsoil=F"
521     intheat=0.
522     call getin_p("intheat",intheat)
523     if (is_master) write(*,*)trim(rname)//": intheat = ",intheat
524
525
526     if (is_master) write(*,*)trim(rname)//": Test physics timescale in 1D?"
527     testradtimes=.false.
528     call getin_p("testradtimes",testradtimes)
529     if (is_master) write(*,*)trim(rname)//": testradtimes = ",testradtimes
530
531! Test of incompatibility:
532! if testradtimes used, we must be in 1D
533     if (testradtimes.and.(ngrid.gt.1)) then
534       call abort_physic(rname,'testradtimes can only be used in 1D!',1)
535     endif
536
537     if (is_master) write(*,*)trim(rname)//": Default planetary temperature?"
538     tplanet=215.0
539     call getin_p("tplanet",tplanet)
540     if (is_master) write(*,*)trim(rname)//": tplanet = ",tplanet
541
542     if (is_master) write(*,*)trim(rname)//": Which star?"
543     startype=1 ! default value = Sol
544     call getin_p("startype",startype)
545     if (is_master) write(*,*)trim(rname)//": startype = ",startype
546
547     if (is_master) write(*,*)trim(rname)//": Value of stellar flux at 1 AU?"
548     Fat1AU=1356.0 ! default value = Sol today
549     call getin_p("Fat1AU",Fat1AU)
550     if (is_master) write(*,*)trim(rname)//": Fat1AU = ",Fat1AU
551
552
553! TRACERS:
554
555     if (is_master) write(*,*)trim(rname)//": fixed radii for Cloud particles?"
556     radfixed=.false. ! default value
557     call getin_p("radfixed",radfixed)
558     if (is_master) write(*,*)trim(rname)//": radfixed = ",radfixed
559
560     if(kastprof)then
561        radfixed=.true.
562     endif
563
564      if (is_master) write(*,*)trim(rname)//&
565       "Number of radiatively active aerosols:"
566     naerkind=0 ! default value
567     call getin_p("naerkind",naerkind)
568     if (is_master) write(*,*)trim(rname)//": naerkind = ",naerkind
569
570
571!***************************************************************
572         !! TRACERS options
573
574     if (is_master)write(*,*)trim(rname)//&
575      "call N2 condensation ?"
576     n2cond=.true. ! default value
577     call getin_p("n2cond",n2cond)
578     if (is_master)write(*,*)trim(rname)//&
579      " n2cond = ",n2cond
580
581     if (is_master)write(*,*)trim(rname)//&
582      "call N2 cloud condensation ?"
583     condensn2=.false. ! default value
584     call getin_p("condensn2",condensn2)
585     if (is_master)write(*,*)trim(rname)//&
586      "condensn2 = ",condensn2
587
588     if (is_master)write(*,*)trim(rname)//&
589      "call no N2 frost formation?"
590     no_n2frost=.false. ! default value
591     call getin_p("no_n2frost",no_n2frost)
592     if (is_master)write(*,*)trim(rname)//&
593      "no_n2frost = ",no_n2frost
594
595     if (is_master)write(*,*)trim(rname)//&
596      "N2 condensation subtimestep?"
597     nbsub=20 ! default value
598     call getin_p("nbsub",nbsub)
599     if (is_master)write(*,*)trim(rname)//&
600      " nbsub = ",nbsub
601
602     if (is_master)write(*,*)trim(rname)//&
603      "Gravitationnal sedimentation ?"
604     sedimentation=.true. ! default value
605     call getin_p("sedimentation",sedimentation)
606     if (is_master)write(*,*)trim(rname)//&
607      " sedimentation = ",sedimentation
608
609     if (is_master)write(*,*)trim(rname)//&
610      "Compute glacial flow ?"
611     glaflow=.false. ! default value
612     call getin_p("glaflow",glaflow)
613     if (is_master)write(*,*)trim(rname)//&
614      " glaflow = ",glaflow
615
616     if (is_master)write(*,*)trim(rname)//&
617      "Compute methane cycle ?"
618     methane=.false. ! default value
619     call getin_p("methane",methane)
620     if (is_master)write(*,*)trim(rname)//&
621      " methane = ",methane
622     condmetsurf=.true. ! default value
623     call getin_p("condmetsurf",condmetsurf)
624     if (is_master)write(*,*)trim(rname)//&
625      " condmetsurf = ",condmetsurf
626     if (is_master)write(*,*)trim(rname)//&
627      "call no CH4 frost formation?"
628     no_ch4frost=.false. ! default value
629     call getin_p("no_ch4frost",no_ch4frost)
630     if (is_master)write(*,*)trim(rname)//&
631      "no_ch4frost = ",no_ch4frost
632
633     if (is_master)write(*,*)trim(rname)//&
634      "Compute methane clouds ?"
635     metcloud=.false. ! default value
636     call getin_p("metcloud",metcloud)
637     if (is_master)write(*,*)trim(rname)//&
638      " metcloud = ",metcloud
639
640     if (is_master)write(*,*)trim(rname)//&
641      "Compute CO cycle ?"
642     carbox=.false. ! default value
643     call getin_p("carbox",carbox)
644     if (is_master)write(*,*)trim(rname)//&
645      " carbox = ",carbox
646     condcosurf=.true. ! default value
647     call getin_p("condcosurf",condcosurf)
648     if (is_master)write(*,*)trim(rname)//&
649      " condcosurf = ",condcosurf
650
651     if (is_master)write(*,*)trim(rname)//&
652      "Compute CO clouds ?"
653     monoxcloud=.false. ! default value
654     call getin_p("monoxcloud",monoxcloud)
655     if (is_master)write(*,*)trim(rname)//&
656      " monoxcloud = ",monoxcloud
657
658     if (is_master)write(*,*)trim(rname)//&
659     "atmospheric redistribution (s):"
660     tau_n2=1. ! default value
661     call getin_p("tau_n2",tau_n2)
662     if (is_master)write(*,*)trim(rname)//&
663     " tau_n2 = ",tau_n2
664     tau_ch4=1.E7 ! default value
665     call getin_p("tau_ch4",tau_ch4)
666     if (is_master)write(*,*)trim(rname)//&
667     " tau_ch4 = ",tau_ch4
668     tau_co=1. ! default value
669     call getin_p("tau_co",tau_co)
670     if (is_master)write(*,*)trim(rname)//&
671     " tau_co = ",tau_co
672     tau_prechaze=1. ! default value
673     call getin_p("tau_prechaze",tau_prechaze)
674     if (is_master)write(*,*)trim(rname)//&
675     " tau_prechaze = ",tau_prechaze
676
677     if (is_master)write(*,*)trim(rname)//&
678      "Day fraction for limited cold trap in SP?"
679     dayfrac=0. ! default value
680     call getin_p("dayfrac",dayfrac)
681     if (is_master)write(*,*)trim(rname)//&
682      " dayfrac = ",dayfrac
683     thresh_non2=0. ! default value
684     call getin_p("thresh_non2",thresh_non2)
685     if (is_master)write(*,*)trim(rname)//&
686      " thresh_non2 = ",thresh_non2
687
688    ! use Pluto.old routines
689     if (is_master) write(*,*) trim(rname)//": use vdifc from old Pluto ?"
690     oldplutovdifc=.false. ! default value
691     call getin_p("oldplutovdifc",oldplutovdifc)
692     if (is_master) write(*,*) trim(rname)//": oldplutovdifc = ",oldplutovdifc
693
694     if (is_master) write(*,*) trim(rname)//&
695       ": call pluto.old correlated-k radiative transfer ?"
696     oldplutocorrk=.false. ! default value
697     call getin_p("oldplutocorrk",oldplutocorrk)
698     if (is_master) write(*,*) trim(rname)//": oldplutocorrk = ",oldplutocorrk
699
700     if (is_master) write(*,*) trim(rname)//&
701       ": call pluto.old sedimentation ?"
702     oldplutosedim=.false. ! default value
703     call getin_p("oldplutosedim",oldplutosedim)
704     if (is_master) write(*,*) trim(rname)//": oldplutosedim = ",oldplutosedim
705
706!***************************************************************
707     !! Haze options
708
709     ! Microphysical moment model
710     ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
711     if (is_master) write(*,*) "Run with or without microphysics?"
712     callmufi=.false. ! default value
713     call getin_p("callmufi",callmufi)
714     if (is_master) write(*,*)" callmufi = ",callmufi
715
716     ! sanity check
717     if (callmufi.and.(.not.tracer)) then
718       print*,"You are running microphysics without tracer"
719       print*,"Please start again with tracer =.true."
720       stop
721     endif
722
723     if (is_master) write(*,*) "Path to microphysical config file?"
724     config_mufi='datagcm/microphysics/config.cfg' ! default value
725     call getin_p("config_mufi",config_mufi)
726     if (is_master) write(*,*)" config_mufi = ",config_mufi
727
728     if (is_master) write(*,*) "Spherical aerosol optical properties datafile"
729     aersprop_file="optprop_rannou_r2-200nm_nu003.dat"  ! default file
730     call getin_p("aersprop_file",aersprop_file)
731     if (is_master) write(*,*) trim(rname)//" aersprop_file = ",trim(aersprop_file)
732
733     if (is_master) write(*,*) "Fractal aerosol optical properties datafile"
734     aerfprop_file="optprop_rannou_fractal_r010nm_N1_1e4_d2.dat"  ! default file
735     call getin_p("aerfprop_file",aerfprop_file)
736     if (is_master) write(*,*) trim(rname)//" aerfprop_file = ",trim(aerfprop_file)
737
738     if (is_master) write(*,*) "Use haze production from CH4 photolysis or production rate?"
739     call_haze_prod_pCH4=.false. ! default value
740     call getin_p("call_haze_prod_pCH4",call_haze_prod_pCH4)
741     if (is_master) write(*,*)" call_haze_prod_pCH4 = ",call_haze_prod_pCH4
742
743     if (is_master) write(*,*) "Pressure level of aerosols production (Pa)?"
744     haze_p_prod=1.0e-2 ! default value
745     call getin_p("haze_p_prod",haze_p_prod)
746     if (is_master) write(*,*)" haze_p_prod = ",haze_p_prod
747
748     if (is_master) write(*,*) "Aerosol production rate (kg.m-2.s-1)?"
749     haze_tx_prod=9.8e-14 ! default value
750     call getin_p("haze_tx_prod",haze_tx_prod)
751     if (is_master) write(*,*)" haze_tx_prod = ",haze_tx_prod
752
753     if (is_master) write(*,*) "Equivalent radius production (m)?"
754     haze_rc_prod=1.0e-9 ! default value
755     call getin_p("haze_rc_prod",haze_rc_prod)
756     if (is_master) write(*,*)" haze_rc_prod = ",haze_rc_prod
757
758     if (is_master) write(*,*) "Monomer radius (m)?"
759     haze_rm=1.0e-8 ! default value
760     call getin_p("haze_rm",haze_rm)
761     if (is_master) write(*,*)" haze_rm = ",haze_rm
762
763     if (is_master) write(*,*) "Aerosol's fractal dimension?"
764     haze_df=2.0 ! default value
765     call getin_p("haze_df",haze_df)
766     if (is_master) write(*,*)" haze_df = ",haze_df
767
768     if (is_master) write(*,*) "Aerosol density (kg.m-3)?"
769     haze_rho=800.0 ! default value
770     call getin_p("haze_rho",haze_rho)
771     if (is_master) write(*,*)" haze_rho = ",haze_rho
772
773     if (is_master) write(*,*) "Radius of air molecule (m)?"
774     air_rad=1.75e-10 ! default value
775     call getin_p("air_rad",air_rad)
776     if (is_master) write(*,*)" air_rad = ",air_rad
777
778     ! Conversion precursor to haze
779     ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
780     if (is_master) write(*,*) "Time constant of conversion in aerosol (s)?"
781     tcon_ch4=1.e7 ! default value
782     call getin_p("tcon_ch4",tcon_ch4)
783     if (is_master) write(*,*)" tcon_ch4 = ",tcon_ch4
784
785     if (is_master) write(*,*) "Ratio CH4 photolysis?"
786     k_ch4=1. ! default value
787     call getin_p("k_ch4",k_ch4)
788     if (is_master) write(*,*)" k_ch4 = ",k_ch4
789
790     if (is_master) write(*,*) "Nitrogen contribution ratio N/C?"
791     ncratio_ch4=0.5 ! default value
792     call getin_p("ncratio_ch4",ncratio_ch4)
793     if (is_master) write(*,*)" ncratio_ch4 = ",ncratio_ch4
794
795     ! Variables for aerosol absorption
796     ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
797     if (is_master) write(*,*) "Sph. aer. absorption correction in VI?"
798     Fabs_aers_VI=1. ! default value
799     call getin_p("Fabs_aers_VI",Fabs_aers_VI)
800     if (is_master) write(*,*)" Fabs_aers_VI = ",Fabs_aers_VI
801
802     if (is_master) write(*,*) "Fra. aer. absorption correction in VI?"
803     Fabs_aerf_VI=1. ! default value
804     call getin_p("Fabs_aerf_VI",Fabs_aerf_VI)
805     if (is_master) write(*,*)" Fabs_aerf_VI = ",Fabs_aerf_VI
806
807     if (is_master) write(*,*) "Sph. aer. absorption correction in IR?"
808     Fabs_aers_IR=1. ! default value
809     call getin_p("Fabs_aers_IR",Fabs_aers_IR)
810     if (is_master) write(*,*)" Fabs_aers_IR = ",Fabs_aers_IR
811
812     if (is_master) write(*,*) "Fra. aer. absorption correction in IR?"
813     Fabs_aerf_IR=1. ! default value
814     call getin_p("Fabs_aerf_IR",Fabs_aerf_IR)
815     if (is_master) write(*,*)" Fabs_aerf_IR = ",Fabs_aerf_IR
816
817     ! Pluto haze model
818     ! ~~~~~~~~~~~~~~~~
819     if (is_master)write(*,*)trim(rname)//&
820     "Production of haze ?"
821     haze=.false. ! default value
822     call getin_p("haze",haze)
823     if (is_master)write(*,*)trim(rname)//&
824     " haze = ",haze
825     hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed
826     call getin_p("hazeconservch4",hazeconservch4)
827     if (is_master)write(*,*)trim(rname)//&
828     "hazconservch4 = ",hazeconservch4
829     if (is_master) write(*,*) "Production of haze with CH4 fix profile?"
830     haze_ch4proffix=.false. ! default value
831     call getin_p("haze_ch4proffix",haze_ch4proffix)
832     if (is_master) write(*,*)" haze_ch4proffix = ",haze_ch4proffix
833     if (is_master)write(*,*)trim(rname)//&
834     "Production of haze (fast version) ?"
835     fasthaze=.false. ! default value
836     call getin_p("fasthaze",fasthaze)
837     if (is_master)write(*,*)trim(rname)//&
838     "fasthaze = ",fasthaze
839
840     if (is_master)write(*,*)trim(rname)//&
841     "Add source of haze ?"
842     source_haze=.false. ! default value
843     call getin_p("source_haze",source_haze)
844     if (is_master)write(*,*)trim(rname)//&
845     " source_haze = ",source_haze
846     mode_hs=0 ! mode haze source default value
847     call getin_p("mode_hs",mode_hs)
848     if (is_master)write(*,*)trim(rname)//&
849     " mode_hs",mode_hs
850     kfix=1 ! default value
851     call getin_p("kfix",kfix)
852     if (is_master)write(*,*)trim(rname)//&
853     "levels kfix",kfix
854     fracsource=0.1 ! default value
855     call getin_p("fracsource",fracsource)
856     if (is_master)write(*,*)trim(rname)//&
857     " fracsource",fracsource
858     latsource=30. ! default value
859     call getin_p("latsource",latsource)
860     if (is_master)write(*,*)trim(rname)//&
861     " latsource",latsource
862     lonsource=180. ! default value
863     call getin_p("lonsource",lonsource)
864     if (is_master)write(*,*)trim(rname)//&
865     " lonsource",lonsource
866
867     if (is_master)write(*,*)trim(rname)//&
868     "Radiatively active haze ?"
869     optichaze=.false. ! default value
870     call getin_p("optichaze",optichaze)
871     if (is_master)write(*,*)trim(rname)//&
872     "optichaze = ",optichaze
873
874     aerohaze=.false. ! default value
875     call getin_p("aerohaze",aerohaze)
876     if (aerohaze) then
877      if (is_master) write(*,*)trim(rname)//": aerohaze is deprecated.",&
878      "it is now called optichaze=.true."
879      call abort_physic(rname,"aerohaze is deprecated. It is now called optichaze",1)
880     endif
881
882     if (is_master)write(*,*)trim(rname)//&
883     "Haze monomer radius ?"
884     rad_haze=10.e-9 ! default value
885     call getin_p("rad_haze",rad_haze)
886     if (is_master)write(*,*)trim(rname)//&
887     "rad_haze = ",rad_haze
888
889     if (is_master)write(*,*)trim(rname)//&
890     "fractal particle ?"
891     fractal=.false. ! default value
892     call getin_p("fractal",fractal)
893     if (is_master)write(*,*)trim(rname)//&
894     "fractal = ",fractal
895     nb_monomer=10 ! default value
896     call getin_p("nb_monomer",nb_monomer)
897     if (is_master)write(*,*)trim(rname)//&
898     "nb_monomer = ",nb_monomer
899
900     if (is_master)write(*,*)trim(rname)//&
901     "fixed gaseous CH4 mixing ratio for rad transfer ?"
902    ch4fix=.false. ! default value
903    call getin_p("ch4fix",ch4fix)
904    if (is_master)write(*,*)trim(rname)//&
905     " ch4fix = ",ch4fix
906    if (is_master)write(*,*)trim(rname)//&
907     "fixed gaseous CH4 mixing ratio for rad transfer ?"
908    vmrch4fix=0.5 ! default value
909    call getin_p("vmrch4fix",vmrch4fix)
910    if (is_master)write(*,*)trim(rname)//&
911     " vmrch4fix = ",vmrch4fix
912    vmrch4_proffix=.false. ! default value
913    call getin_p("vmrch4_proffix",vmrch4_proffix)
914    if (is_master)write(*,*)trim(rname)//&
915     " vmrch4_proffix = ",vmrch4_proffix
916
917    if (is_master)write(*,*)trim(rname)//&
918     "call specific cooling for radiative transfer ?"
919    cooling=.false.  ! default value
920    call getin_p("cooling",cooling)
921    if (is_master)write(*,*)trim(rname)//&
922     " cooling = ",cooling
923
924    if (is_master)write(*,*)trim(rname)//&
925     "NLTE correction for methane heating rates?"
926    nlte=.false.  ! default value
927    call getin_p("nlte",nlte)
928    if (is_master)write(*,*)trim(rname)//&
929     " nlte = ",nlte
930    strobel=.true.  ! default value
931    call getin_p("strobel",strobel)
932    if (is_master)write(*,*)trim(rname)//&
933     " strobel = ",strobel
934
935     if (is_master)write(*,*)trim(rname)//&
936     "fixed radius profile from txt file ?"
937     haze_radproffix=.false. ! default value
938     call getin_p("haze_radproffix",haze_radproffix)
939     if (is_master)write(*,*)trim(rname)//&
940     "haze_radproffix = ",haze_radproffix
941     if (is_master)write(*,*)trim(rname)//&
942     "fixed MMR profile from txt file ?"
943     haze_proffix=.false. ! default value
944     call getin_p("haze_proffix",haze_proffix)
945     if (is_master)write(*,*)trim(rname)//&
946     "haze_proffix = ",haze_proffix
947
948     if (is_master)write(*,*)trim(rname)//&
949     "Number mix ratio of haze particles for co clouds:"
950     Nmix_co=100000. ! default value
951     call getin_p("Nmix_co",Nmix_co)
952     if (is_master)write(*,*)trim(rname)//&
953     " Nmix_co = ",Nmix_co
954
955     if (is_master)write(*,*)trim(rname)//&
956     "Number mix ratio of haze particles for ch4 clouds:"
957     Nmix_ch4=100000. ! default value
958     call getin_p("Nmix_ch4",Nmix_ch4)
959     if (is_master)write(*,*)trim(rname)//&
960     " Nmix_ch4 = ",Nmix_ch4
961
962!***************************************************************
963     !! Surface properties
964
965!*********** N2 *********************************
966
967     if (is_master)write(*,*)trim(rname)//&
968      "Mode for changing N2 albedo"
969     mode_n2=0 ! default value
970     call getin_p("mode_n2",mode_n2)
971     if (is_master)write(*,*)trim(rname)//&
972      " mode_n2 = ",mode_n2
973     thres_n2ice=1. ! default value
974     call getin_p("thres_n2ice",thres_n2ice)
975     if (is_master)write(*,*)trim(rname)//&
976      " thres_n2ice = ",thres_n2ice
977
978     if (is_master)write(*,*)trim(rname)//&
979      "Diff of N2 albedo with thickness"
980     deltab=0. ! default value
981     call getin_p("deltab",deltab)
982     if (is_master)write(*,*)trim(rname)//&
983      " deltab = ",deltab
984
985     if (is_master)write(*,*)trim(rname)//&
986      "albedo N2 beta "
987     alb_n2b=0.7 ! default value
988     call getin_p("alb_n2b",alb_n2b)
989     if (is_master)write(*,*)trim(rname)//&
990      " alb_n2b = ",alb_n2b
991
992     if (is_master)write(*,*)trim(rname)//&
993      "albedo N2 alpha "
994     alb_n2a=0.7 ! default value
995     call getin_p("alb_n2a",alb_n2a)
996     if (is_master)write(*,*)trim(rname)//&
997      " alb_n2a = ",alb_n2a
998
999     if (is_master)write(*,*)trim(rname)//&
1000      "emis N2 beta "
1001     emis_n2b=0.7 ! default value
1002     call getin_p("emis_n2b",emis_n2b)
1003     if (is_master)write(*,*)trim(rname)//&
1004      " emis_n2b = ",emis_n2b
1005
1006     if (is_master)write(*,*)trim(rname)//&
1007      "emis N2 alpha "
1008     emis_n2a=0.7 ! default value
1009     call getin_p("emis_n2a",emis_n2a)
1010     if (is_master)write(*,*)trim(rname)//&
1011      " emis_n2a = ",emis_n2a
1012
1013!*********** CH4 *********************************
1014
1015     if (is_master)write(*,*)trim(rname)//&
1016      "Mode for changing CH4 albedo"
1017     mode_ch4=0 ! default value
1018     call getin_p("mode_ch4",mode_ch4)
1019     if (is_master)write(*,*)trim(rname)//&
1020      " mode_ch4 = ",mode_ch4
1021     feedback_met=0 ! default value
1022     call getin_p("feedback_met",feedback_met)
1023     if (is_master)write(*,*)trim(rname)//&
1024      " feedback_met = ",feedback_met
1025     thres_ch4ice=1. ! default value, mm
1026     call getin_p("thres_ch4ice",thres_ch4ice)
1027     if (is_master)write(*,*)trim(rname)//&
1028      " thres_ch4ice = ",thres_ch4ice
1029     fdch4_latn=-1.
1030     fdch4_lats=0.
1031     fdch4_lone=0.
1032     fdch4_lonw=-1.
1033     fdch4_depalb=0.5
1034     fdch4_finalb=0.95
1035     fdch4_maxalb=0.99
1036     fdch4_ampl=200.
1037     fdch4_maxice=100.
1038     call getin_p("fdch4_latn",fdch4_latn)
1039     call getin_p("fdch4_lats",fdch4_lats)
1040     call getin_p("fdch4_lone",fdch4_lone)
1041     call getin_p("fdch4_lonw",fdch4_lonw)
1042     call getin_p("fdch4_depalb",fdch4_depalb)
1043     call getin_p("fdch4_finalb",fdch4_finalb)
1044     call getin_p("fdch4_maxalb",fdch4_maxalb)
1045     call getin_p("fdch4_maxice",fdch4_maxice)
1046     call getin_p("fdch4_ampl",fdch4_ampl)
1047     if (is_master)write(*,*)trim(rname)//&
1048      " Values for albedo feedback = ",fdch4_latn,&
1049     fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,&
1050     fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl
1051
1052     if (is_master)write(*,*)trim(rname)//&
1053      "Latitude for diff albedo methane"
1054     metlateq=25. ! default value
1055     call getin_p("metlateq",metlateq)
1056     if (is_master)write(*,*)trim(rname)//&
1057      " metlateq = ",metlateq
1058
1059     if (is_master)write(*,*)trim(rname)//&
1060      "Ls1 and Ls2 for change of ch4 albedo"
1061     metls1=-1. ! default value
1062     metls2=-2. ! default value
1063     call getin_p("metls1",metls1)
1064     call getin_p("metls2",metls2)
1065
1066     if (is_master)write(*,*)trim(rname)//&
1067      "albedo CH4 "
1068     alb_ch4=0.5 ! default value
1069     call getin_p("alb_ch4",alb_ch4)
1070     if (is_master)write(*,*)trim(rname)//&
1071      " alb_ch4 = ",alb_ch4
1072
1073     if (is_master)write(*,*)trim(rname)//&
1074      "albedo equatorial CH4 "
1075     alb_ch4_eq=alb_ch4 ! default value
1076     call getin_p("alb_ch4_eq",alb_ch4_eq)
1077     if (is_master)write(*,*)trim(rname)//&
1078      " alb_ch4_eq = ",alb_ch4_eq
1079
1080     if (is_master)write(*,*)trim(rname)//&
1081      "albedo south hemis CH4 "
1082     alb_ch4_s=alb_ch4 ! default value
1083     call getin_p("alb_ch4_s",alb_ch4_s)
1084     if (is_master)write(*,*)trim(rname)//&
1085      " alb_ch4_s = ",alb_ch4_s
1086
1087     if (is_master)write(*,*)trim(rname)//&
1088      "emis CH4 "
1089     emis_ch4=0.5 ! default value
1090     call getin_p("emis_ch4",emis_ch4)
1091     if (is_master)write(*,*)trim(rname)//&
1092      " emis_ch4 = ",emis_ch4
1093
1094     if (is_master)write(*,*)trim(rname)//&
1095      "CH4 lag for n2 sublimation limitation"
1096     ch4lag=.false. ! default value
1097     latlag=-90. ! default value
1098     vmrlag=1. ! default value
1099     call getin_p("ch4lag",ch4lag)
1100     call getin_p("latlag",latlag)
1101     call getin_p("vmrlag",vmrlag)
1102     if (is_master)write(*,*)trim(rname)//&
1103      " ch4lag = ",ch4lag
1104     if (is_master)write(*,*)trim(rname)//&
1105      " latlag = ",latlag
1106     if (is_master)write(*,*)trim(rname)//&
1107      " vmrlag = ",vmrlag
1108
1109     if (is_master)write(*,*)trim(rname)//&
1110      "Max temperature for surface ?"
1111     tsurfmax=.false. ! default value
1112     albmin_ch4=0.3 ! default value
1113     call getin_p("tsurfmax",tsurfmax)
1114     call getin_p("albmin_ch4",albmin_ch4)
1115     if (is_master)write(*,*)trim(rname)//&
1116      " tsurfmax = ",tsurfmax
1117     if (is_master)write(*,*)trim(rname)//&
1118      " albmin_ch4 = ",albmin_ch4
1119
1120
1121     if (is_master)write(*,*)trim(rname)//&
1122     "fixed gaseous CH4 mixing ratio for rad transfer ?"
1123     ch4fix=.false. ! default value
1124     call getin_p("ch4fix",ch4fix)
1125     if (is_master)write(*,*)trim(rname)//&
1126       " ch4fix = ",ch4fix
1127     if (is_master)write(*,*)trim(rname)//&
1128       "fixed gaseous CH4 mixing ratio for rad transfer ?"
1129     vmrch4fix=0.5 ! default value
1130     call getin_p("vmrch4fix",vmrch4fix)
1131     if (is_master)write(*,*)trim(rname)//&
1132        " vmrch4fix = ",vmrch4fix
1133     vmrch4_proffix=.false. ! default value
1134     call getin_p("vmrch4_proffix",vmrch4_proffix)
1135     if (is_master)write(*,*)trim(rname)//&
1136        " vmrch4_proffix = ",vmrch4_proffix
1137
1138
1139!*********** CO *********************************
1140
1141     if (is_master)write(*,*)trim(rname)//&
1142      "albedo CO "
1143     alb_co=0.4 ! default value
1144     call getin_p("alb_co",alb_co)
1145     if (is_master)write(*,*)trim(rname)//&
1146      " alb_co = ",alb_co
1147     thres_coice=1. ! default value, mm
1148     call getin_p("thres_coice",thres_coice)
1149     if (is_master)write(*,*)trim(rname)//&
1150      " thres_coice = ",thres_coice
1151
1152     if (is_master)write(*,*)trim(rname)//&
1153      "emis CO "
1154     emis_co=0.4 ! default value
1155     call getin_p("emis_co",emis_co)
1156     if (is_master)write(*,*)trim(rname)//&
1157      " emis_co = ",emis_co
1158
1159!*********** THOLINS *********************************
1160     if (is_master)write(*,*)trim(rname)//&
1161      "Mode for changing tholins albedo/emis"
1162     mode_tholins=0 ! default value
1163     call getin_p("mode_tholins",mode_tholins)
1164     if (is_master)write(*,*)trim(rname)//&
1165      " mode_tholins = ",mode_tholins
1166
1167     if (is_master)write(*,*)trim(rname)//&
1168      "albedo tho "
1169     alb_tho=0.1 ! default value
1170     call getin_p("alb_tho",alb_tho)
1171     if (is_master)write(*,*)trim(rname)//&
1172      " alb_tho = ",alb_tho
1173
1174     if (is_master)write(*,*)trim(rname)//&
1175      "albedo tho eq"
1176     alb_tho_eq=0.1 ! default value
1177     call getin_p("alb_tho_eq",alb_tho_eq)
1178     if (is_master)write(*,*)trim(rname)//&
1179      " alb_tho_eq = ",alb_tho_eq
1180
1181     if (is_master)write(*,*)trim(rname)//&
1182      "emis tholins "
1183     emis_tho=1. ! default value
1184     call getin_p("emis_tho",emis_tho)
1185     if (is_master)write(*,*)trim(rname)//&
1186      " emis_tho = ",emis_tho
1187
1188     if (is_master)write(*,*)trim(rname)//&
1189      "emis tholins eq"
1190     emis_tho_eq=1. ! default value
1191     call getin_p("emis_tho_eq",emis_tho_eq)
1192     if (is_master)write(*,*)trim(rname)//&
1193      " emis_tho_eq = ",emis_tho_eq
1194
1195     if (is_master)write(*,*)trim(rname)//&
1196      "Latitude for diff albedo tholins"
1197     tholateq=25. ! default value
1198     call getin_p("tholateq",tholateq)
1199     if (is_master)write(*,*)trim(rname)//&
1200      " tholateq = ",tholateq
1201     tholatn=-1.
1202     tholats=0.
1203     tholone=0.
1204     tholonw=-1.
1205     alb_tho_spe=0.1 ! default value
1206     emis_tho_spe=1. ! default value
1207     call getin_p("  tholatn",tholatn)
1208     call getin_p("  tholats",tholats)
1209     call getin_p("  tholone",tholone)
1210     call getin_p("  tholonw",tholonw)
1211     if (is_master)write(*,*)trim(rname)//&
1212      " Values for special tholins albedo = ",tholatn,&
1213        tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe
1214
1215     if (is_master)write(*,*)trim(rname)//&
1216      "Specific albedo"
1217     spelon1=-180. ! default value
1218     spelon2=180. ! default value
1219     spelat1=-90. ! default value
1220     spelat2=90. ! default value
1221     specalb=.false. ! default value
1222     if (is_master)write(*,*)trim(rname)//&
1223      "albedo/emis spe "
1224     albspe=0.1 ! default value
1225     emispe=1. ! default value
1226     call getin_p("spelon1",spelon1)
1227     call getin_p("spelon2",spelon2)
1228     call getin_p("spelat1",spelat1)
1229     call getin_p("spelat2",spelat2)
1230     call getin_p("specalb",specalb)
1231     call getin_p("albspe",albspe)
1232     call getin_p("emispe",emispe)
1233
1234     if (is_master)write(*,*)trim(rname)//&
1235      " specific = ",specalb
1236
1237!********************** TI *****************************
1238
1239     if (is_master)write(*,*)trim(rname)//&
1240      "Change TI with time"
1241     changeti=.false. ! default value
1242     call getin_p("changeti",changeti)
1243     if (is_master)write(*,*)trim(rname)//&
1244      " changeti = ",changeti
1245     changetid=.false. ! default value for diurnal TI
1246     call getin_p("changetid",changetid)
1247     if (is_master)write(*,*)trim(rname)//&
1248      " changetid = ",changetid
1249
1250     if (is_master)write(*,*)trim(rname)//&
1251      "IT N2 "
1252     ITN2=800. ! default value
1253     call getin_p("ITN2",ITN2)
1254     if (is_master)write(*,*)trim(rname)//&
1255      " ITN2 = ",ITN2
1256     ITN2d=20. ! default value
1257     call getin_p("ITN2d",ITN2d)
1258     if (is_master)write(*,*)trim(rname)//&
1259      " ITN2d = ",ITN2d
1260
1261     if (is_master)write(*,*)trim(rname)//&
1262      "IT CH4"
1263     ITCH4=800. ! default value
1264     call getin_p("ITCH4",ITCH4)
1265     if (is_master)write(*,*)trim(rname)//&
1266      " ITCH4 = ",ITCH4
1267     ITCH4d=20. ! default value
1268     call getin_p("ITCH4d",ITCH4d)
1269     if (is_master)write(*,*)trim(rname)//&
1270      " ITCH4d = ",ITCH4d
1271
1272     if (is_master)write(*,*)trim(rname)//&
1273      "IT H2O"
1274     ITH2O=800. ! default value
1275     call getin_p("ITH2O",ITH2O)
1276     if (is_master)write(*,*)trim(rname)//&
1277      " ITH2O = ",ITH2O
1278     ITH2Od=20. ! default value
1279     call getin_p("ITH2Od",ITH2Od)
1280     if (is_master)write(*,*)trim(rname)//&
1281      " ITH2Od = ",ITH2Od
1282
1283!************** COOLING ***************
1284
1285     alpha_top=5e-11 ! default value
1286     call getin_p("alpha_top",alpha_top)
1287     if (is_master)write(*,*)trim(rname)//&
1288      " alpha_top = ",alpha_top
1289     pref=0.02 ! default value
1290     call getin_p("pref",pref)
1291     if (is_master)write(*,*)trim(rname)//&
1292      " pref = ",pref
1293     deltap=0.1 ! default value
1294     call getin_p("deltap",deltap)
1295     if (is_master)write(*,*)trim(rname)//&
1296      " deltap = ",deltap
1297
1298!=================================
1299
1300     if (is_master) write(*,*)trim(rname)//&
1301       ": Is the variable gas species radiatively active?"
1302     Tstrat=167.0
1303     varactive=.false.
1304     call getin_p("varactive",varactive)
1305     if (is_master) write(*,*)trim(rname)//": varactive = ",varactive
1306
1307     if (is_master) write(*,*)trim(rname)//&
1308       ": Is the variable gas species distribution set?"
1309     varfixed=.false.
1310     call getin_p("varfixed",varfixed)
1311     if (is_master) write(*,*)trim(rname)//": varfixed = ",varfixed
1312
1313! Test of incompatibility:
1314! if varactive, then varfixed should be false
1315     if (varactive.and.varfixed) then
1316       call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
1317     endif
1318
1319     if (is_master) write(*,*)trim(rname)//": Gravitationnal sedimentation ?"
1320     sedimentation=.false. ! default value
1321     call getin_p("sedimentation",sedimentation)
1322     if (is_master) write(*,*)trim(rname)//": sedimentation = ",sedimentation
1323
1324     if (is_master) write(*,*)trim(rname)//": Spectral Dependant albedo ?"
1325     albedo_spectral_mode=.false. ! default value
1326     call getin_p("albedo_spectral_mode",albedo_spectral_mode)
1327     if (is_master) write(*,*)trim(rname)//&
1328     ": albedo_spectral_mode = ",albedo_spectral_mode
1329
1330     if (is_master) then
1331       write(*,*)trim(rname)//": Snow albedo ?"
1332       write(*,*)trim(rname)//": If albedo_spectral_mode=.true., then this "
1333       write(*,*)trim(rname)//": corresponds to the 0.5 microns snow albedo."
1334     endif
1335     albedosnow=0.5         ! default value
1336     call getin_p("albedosnow",albedosnow)
1337     if (is_master) write(*,*)trim(rname)//": albedosnow = ",albedosnow
1338
1339     if (is_master) write(*,*)trim(rname)//": N2 ice albedo ?"
1340     albedon2ice=0.5       ! default value
1341     call getin_p("albedon2ice",albedon2ice)
1342     if (is_master) write(*,*)trim(rname)//": albedon2ice = ",albedon2ice
1343
1344     if (is_master) write(*,*)trim(rname)//": Maximum ice thickness ?"
1345     maxicethick=2.0         ! default value
1346     call getin_p("maxicethick",maxicethick)
1347     if (is_master) write(*,*)trim(rname)//": maxicethick = ",maxicethick
1348
1349     if (is_master) write(*,*)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1350     kmixmin=1.0e-2         ! default value
1351     call getin_p("kmixmin",kmixmin)
1352     if (is_master) write(*,*)trim(rname)//": kmixmin = ",kmixmin
1353
1354     if (is_master) write(*,*)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
1355     if (is_master) write(*,*)'Predefined Mg from dynamics is ',mugaz,'amu'
1356
1357     force_cpp=.false. ! default value
1358     call getin_p("force_cpp",force_cpp)
1359     if (force_cpp) then
1360      if (is_master) write(*,*)trim(rname)//": force_cpp = ",force_cpp
1361      if (is_master) write(*,*)trim(rname)//": force_cpp is deprecated.",&
1362      "Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
1363      call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
1364     endif
1365
1366     if (is_master) write(*,*)trim(rname)//&
1367     ": where do you want your cpp/mugaz value to come from?",&
1368     "=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
1369     cpp_mugaz_mode = 0 ! default value
1370     call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
1371     if (is_master) write(*,*)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
1372
1373
1374! Test of incompatibility:
1375
1376     if ((.not.tracer).and.(haze)) then
1377       call abort_physic(rname, 'if haze are on, tracers must be on!', 1)
1378     endif
1379     if ((callmufi).and.(haze)) then
1380       call abort_physic(rname, 'if haze are on, microphysics should be deactivated!', 1)
1381     endif
1382     if ((haze).and.(naerkind.gt.1)) then
1383      call abort_physic(rname, 'if haze are on, naerkind must be equal to 1!', 1)
1384     endif
1385     if ((callmufi).and..not.(naerkind.gt.1)) then
1386      call abort_physic(rname, 'if microphysics is on, naerkind must be > 1!', 1)
1387     endif
1388     if (.not.(callmufi.or.haze).and.(optichaze)) then
1389      call abort_physic(rname, 'if microphysics and haze are off, optichaze must be deactivated!', 1)
1390     endif
1391     if ((callmufi.and.call_haze_prod_pCH4).and..not.(methane)) then
1392      call abort_physic(rname, 'if haze production from CH4 photolysis is on, methane must be activated!', 1)
1393     endif
1394     if (haze_proffix.and.sedimentation) then
1395         call abort_physic(rname, 'if haze profile is set, sedimentation must be deactivated', 1)
1396     endif
1397     if (callmolvis.and..not.callconduct) then
1398         call abort_physic(rname, 'if callmolvis is set, callconduct must be true', 1)
1399     endif
1400     if (glaflow.and..not.fast) then
1401         call abort_physic(rname, 'if glaflow is set, fast must be true', 1)
1402     endif
1403     if (paleo.and..not.fast) then
1404         call abort_physic(rname, 'if paleo is set, fast must be true', 1)
1405     endif
1406     if ((haze_proffix.or.haze_radproffix).and..not.optichaze) then
1407         call abort_physic(rname, 'for now, haze/rad_proffix only works with optichaze=T', 1)
1408     endif
1409     ! if (carbox.and.condcosurf.and.no_n2frost) then
1410     !   call abort_physic(rname, "CO surface condensation and no_n2frost are both active which may not be relevant", 1)
1411     ! end if
1412
1413     if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
1414        write(*,*)'    !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
1415        write(*,*)'    in 3D can result in very pathological behavior. You have been warned !!!!!'
1416      else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
1417        ! for this specific 1D error we will remove run.def before aborting JL22
1418        call system("rm -rf run.def")
1419        call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
1420     endif
1421
1422     if (cpp_mugaz_mode == 1) then
1423       mugaz = -99999.
1424       if (is_master) write(*,*)trim(rname)//&
1425         ": MEAN MOLECULAR MASS in g mol-1 ?"
1426       call getin_p("mugaz",mugaz)
1427       IF (mugaz.eq.-99999.) THEN
1428         call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
1429       ENDIF
1430       cpp = -99999.
1431       if (is_master) write(*,*)trim(rname)//&
1432         ": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1433       call getin_p("cpp",cpp)
1434       IF (cpp.eq.-99999.) THEN
1435           PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
1436           STOP
1437       ENDIF
1438       if (is_master) write(*,*)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
1439       if (is_master) write(*,*)'New Mg from callfis is ',mugaz,'amu'
1440
1441     endif ! of if (cpp_mugaz_mode == 1)
1442     call su_gases
1443     call calc_cpp_mugaz
1444
1445     if (is_master) then
1446       PRINT*,'--------------------------------------------'
1447       PRINT*
1448       PRINT*
1449     endif
1450  ELSE
1451     call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1452  ENDIF ! of IF(iscallphys)
1453
1454  if (is_master) then
1455    PRINT*
1456    PRINT*,'inifis: daysec',daysec
1457    PRINT*
1458    PRINT*,'inifis: The radiative transfer is computed:'
1459    PRINT*,'           each ',iradia,' physical time-step'
1460    PRINT*,'        or each ',iradia*dtphys,' seconds'
1461    PRINT*
1462  endif
1463
1464!-----------------------------------------------------------------------
1465!     Some more initialization:
1466!     ------------------------
1467
1468  ! Initializations for comgeomfi_h
1469#ifndef MESOSCALE
1470  totarea=SSUM(ngrid,parea,1)
1471  call planetwide_sumval(parea,totarea_planet)
1472
1473  !! those are defined in comdiurn_h.F90
1474  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1475  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1476  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1477  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1478
1479  DO ig=1,ngrid
1480     sinlat(ig)=sin(plat(ig))
1481     coslat(ig)=cos(plat(ig))
1482     sinlon(ig)=sin(plon(ig))
1483     coslon(ig)=cos(plon(ig))
1484  ENDDO
1485#endif
1486  ! initialize variables in radinc_h
1487  call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1488
1489  ! initialize variables and allocate arrays in radcommon_h
1490  call ini_radcommon_h(naerkind)
1491
1492  ! allocate "comsoil_h" arrays
1493  call ini_comsoil_h(ngrid)
1494
1495  END SUBROUTINE inifis
1496
1497END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.