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

Last change on this file since 3652 was 3652, checked in by tbertrand, 4 months ago

Pluto : for fast (nogcm) option only : adding possibility to run the CH4 cycle without any formation of CH4 frost (sublimation/condensation active over perenne reservoirs only)
TB

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