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

Last change on this file since 1672 was 1672, checked in by jvatant, 8 years ago

Re-plug chemistry of old Titan
+ minor correction on reading Titan's startfiles and haze routine
JVO

File size: 17.7 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 radinc_h, only: ini_radinc_h, naerkind
12  use radcommon_h, only: ini_radcommon_h
13  use datafile_mod, only: datadir
14  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
15  use comgeomfi_h, only: totarea, totarea_planet
16  use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
17  use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
18                              init_time, daysec, dtphys
19  use comcstfi_mod, only: rad, cpp, g, r, rcp, &
20                          mugaz, pi, avocado, kbol
21  use planete_mod, only: nres
22  use planetwide_mod, only: planetwide_sumval
23  use callkeys_mod
24  use mod_phys_lmdz_para, only : is_parallel
25
26!=======================================================================
27!
28!   purpose:
29!   -------
30!
31!   Initialisation for the physical parametrisations of the LMD
32!   Generic Model.
33!
34!   author: Frederic Hourdin 15 / 10 /93
35!   -------
36!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
37!             Ehouarn Millour (oct. 2008) tracers are now identified
38!              by their names and may not be contiguously
39!              stored in the q(:,:,:,:) array
40!             E.M. (june 2009) use getin routine to load parameters
41!
42!
43!   arguments:
44!   ----------
45!
46!   input:
47!   ------
48!
49!    ngrid                 Size of the horizontal grid.
50!                          All internal loops are performed on that grid.
51!    nlayer                Number of vertical layers.
52!    pdayref               Day of reference for the simulation
53!    pday                  Number of days counted from the North. Spring
54!                          equinoxe.
55!
56!=======================================================================
57!
58!-----------------------------------------------------------------------
59!   declarations:
60!   -------------
61  use datafile_mod, only: datadir
62  use ioipsl_getin_p_mod, only: getin_p
63  IMPLICIT NONE
64
65
66
67  REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
68  INTEGER,INTENT(IN) :: nday
69  INTEGER,INTENT(IN) :: ngrid,nlayer,nq
70  REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
71  integer,intent(in) :: day_ini
72  INTEGER ig,ierr
73 
74  EXTERNAL iniorbit,orbite
75  EXTERNAL SSUM
76  REAL SSUM
77 
78  ! initialize constants in comcstfi_mod
79  rad=prad
80  cpp=pcpp
81  g=pg
82  r=pr
83  rcp=r/cpp
84  pi=2.*asin(1.)
85  avocado = 6.02214179e23   ! added by RW
86  kbol = 1.38064852e-23  ! added by JVO for Titan chem
87
88  ! Initialize some "temporal and calendar" related variables
89  CALL init_time(day_ini,pdaysec,nday,ptimestep)
90
91  ! read in some parameters from "run.def" for physics,
92  ! or shared between dynamics and physics.
93  call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
94                                    ! in dynamical steps
95  call getin_p("day_step",day_step) ! number of dynamical steps per day
96  call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
97
98  ! do we read a startphy.nc file? (default: .true.)
99  call getin_p("startphy_file",startphy_file)
100 
101! --------------------------------------------------------------
102!  Reading the "callphys.def" file controlling some key options
103! --------------------------------------------------------------
104     
105  ! check that 'callphys.def' file is around
106  OPEN(99,file='callphys.def',status='old',form='formatted',iostat=ierr)
107  CLOSE(99)
108  IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module
109     
110!!!      IF(ierr.EQ.0) THEN
111  IF(iscallphys) THEN
112     PRINT*
113     PRINT*
114     PRINT*,'--------------------------------------------'
115     PRINT*,' inifis: Parametres pour la physique (callphys.def)'
116     PRINT*,'--------------------------------------------'
117
118     write(*,*) "Directory where external input files are:"
119     ! default 'datadir' is set in "datadir_mod"
120     call getin_p("datadir",datadir) ! default path
121     write(*,*) " datadir = ",trim(datadir)
122
123     write(*,*) "Run with or without tracer transport ?"
124     tracer=.false. ! default value
125     call getin_p("tracer",tracer)
126     write(*,*) " tracer = ",tracer
127
128     write(*,*) "Run with or without atm mass update ", &
129            " due to tracer evaporation/condensation?"
130     mass_redistrib=.false. ! default value
131     call getin_p("mass_redistrib",mass_redistrib)
132     write(*,*) " mass_redistrib = ",mass_redistrib
133
134     write(*,*) "Diurnal cycle ?"
135     write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
136     diurnal=.true. ! default value
137     call getin_p("diurnal",diurnal)
138     write(*,*) " diurnal = ",diurnal
139
140     write(*,*) "Seasonal cycle ?"
141     write(*,*) "(if season=false, Ls stays constant, to value ", &
142         "set in 'start'"
143     season=.true. ! default value
144     call getin_p("season",season)
145     write(*,*) " season = ",season
146
147     write(*,*) "Tidally resonant rotation ?"
148     tlocked=.false. ! default value
149     call getin_p("tlocked",tlocked)
150     write(*,*) "tlocked = ",tlocked
151
152     write(*,*) "Saturn ring shadowing ?"
153     rings_shadow = .false.
154     call getin_p("rings_shadow", rings_shadow)
155     write(*,*) "rings_shadow = ", rings_shadow
156         
157     write(*,*) "Compute latitude-dependent gravity field?"
158     oblate = .false.
159     call getin_p("oblate", oblate)
160     write(*,*) "oblate = ", oblate
161
162     write(*,*) "Flattening of the planet (a-b)/a "
163     flatten = 0.0
164     call getin_p("flatten", flatten)
165     write(*,*) "flatten = ", flatten         
166
167     write(*,*) "Needed if oblate=.true.: J2"
168     J2 = 0.0
169     call getin_p("J2", J2)
170     write(*,*) "J2 = ", J2
171         
172     write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
173     MassPlanet = 0.0
174     call getin_p("MassPlanet", MassPlanet)
175     write(*,*) "MassPlanet = ", MassPlanet         
176
177     write(*,*) "Needed if oblate=.true.: Planet mean radius (m)"
178     Rmean = 0.0
179     call getin_p("Rmean", Rmean)
180     write(*,*) "Rmean = ", Rmean
181         
182! Test of incompatibility:
183! if tlocked, then diurnal should be false
184     if (tlocked.and.diurnal) then
185       print*,'If diurnal=true, we should turn off tlocked.'
186       stop
187     endif
188
189     write(*,*) "Tidal resonance ratio ?"
190     nres=0          ! default value
191     call getin_p("nres",nres)
192     write(*,*) "nres = ",nres
193
194     write(*,*) "Write some extra output to the screen ?"
195     lwrite=.false. ! default value
196     call getin_p("lwrite",lwrite)
197     write(*,*) " lwrite = ",lwrite
198
199     write(*,*) "Save statistics in file stats.nc ?"
200     callstats=.true. ! default value
201     call getin_p("callstats",callstats)
202     write(*,*) " callstats = ",callstats
203
204     write(*,*) "Test energy conservation of model physics ?"
205     enertest=.false. ! default value
206     call getin_p("enertest",enertest)
207     write(*,*) " enertest = ",enertest
208
209     write(*,*) "Check to see if cpp values used match gases.def ?"
210     check_cpp_match=.true. ! default value
211     call getin_p("check_cpp_match",check_cpp_match)
212     write(*,*) " check_cpp_match = ",check_cpp_match
213
214     write(*,*) "call radiative transfer ?"
215     callrad=.true. ! default value
216     call getin_p("callrad",callrad)
217     write(*,*) " callrad = ",callrad
218
219     write(*,*) "call correlated-k radiative transfer ?"
220     corrk=.true. ! default value
221     call getin_p("corrk",corrk)
222     write(*,*) " corrk = ",corrk
223
224     write(*,*) "prohibit calculations outside corrk T grid?"
225     strictboundcorrk=.true. ! default value
226     call getin_p("strictboundcorrk",strictboundcorrk)
227     write(*,*) "strictboundcorrk = ",strictboundcorrk
228
229     write(*,*) "call gaseous absorption in the visible bands?", &
230                    "(matters only if callrad=T)"
231     callgasvis=.false. ! default value
232     call getin_p("callgasvis",callgasvis)
233     write(*,*) " callgasvis = ",callgasvis
234       
235     write(*,*) "call continuum opacities in radiative transfer ?", &
236                    "(matters only if callrad=T)"
237     continuum=.true. ! default value
238     call getin_p("continuum",continuum)
239     write(*,*) " continuum = ",continuum
240 
241     write(*,*) "call turbulent vertical diffusion ?"
242     calldifv=.true. ! default value
243     call getin_p("calldifv",calldifv)
244     write(*,*) " calldifv = ",calldifv
245
246     write(*,*) "use turbdiff instead of vdifc ?"
247     UseTurbDiff=.true. ! default value
248     call getin_p("UseTurbDiff",UseTurbDiff)
249     write(*,*) " UseTurbDiff = ",UseTurbDiff
250
251     write(*,*) "call convective adjustment ?"
252     calladj=.true. ! default value
253     call getin_p("calladj",calladj)
254     write(*,*) " calladj = ",calladj
255
256     write(*,*) "Radiative timescale for Newtonian cooling ?"
257     tau_relax=30. ! default value
258     call getin_p("tau_relax",tau_relax)
259     write(*,*) " tau_relax = ",tau_relax
260     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
261
262     write(*,*)"call thermal conduction in the soil ?"
263     callsoil=.true. ! default value
264     call getin_p("callsoil",callsoil)
265     write(*,*) " callsoil = ",callsoil
266         
267     write(*,*)"Rad transfer is computed every iradia", &
268                   " physical timestep"
269     iradia=1 ! default value
270     call getin_p("iradia",iradia)
271     write(*,*)" iradia = ",iradia
272       
273     write(*,*)"Rayleigh scattering ?"
274     rayleigh=.false.
275     call getin_p("rayleigh",rayleigh)
276     write(*,*)" rayleigh = ",rayleigh
277
278     write(*,*) "Use blackbody for stellar spectrum ?"
279     stelbbody=.false. ! default value
280     call getin_p("stelbbody",stelbbody)
281     write(*,*) " stelbbody = ",stelbbody
282
283     write(*,*) "Stellar blackbody temperature ?"
284     stelTbb=5800.0 ! default value
285     call getin_p("stelTbb",stelTbb)
286     write(*,*) " stelTbb = ",stelTbb
287
288     write(*,*)"Output mean OLR in 1D?"
289     meanOLR=.false.
290     call getin_p("meanOLR",meanOLR)
291     write(*,*)" meanOLR = ",meanOLR
292
293     write(*,*)"Output spectral OLR in 3D?"
294     specOLR=.false.
295     call getin_p("specOLR",specOLR)
296     write(*,*)" specOLR = ",specOLR
297
298     write(*,*)"Uniform absorption in radiative transfer?"
299     graybody=.false.
300     call getin_p("graybody",graybody)
301     write(*,*)" graybody = ",graybody
302
303! Chemistry
304
305     write(*,*) "Run with or without chemistry?"
306     callchim=.false. ! default value
307     call getin_p("callchim",callchim)
308     write(*,*) " callchim = ",callchim
309
310     ! sanity check
311     if (callchim.and.(.not.tracer)) then
312       print*,"You are running chemistry without tracer"
313       print*,"Please start again with tracer =.true."
314       stop
315     endif
316
317     write(*,*)"Chemistry is computed every ichim", &
318                   " physical timestep"
319     ichim=1 ! default value
320     call getin_p("ichim",ichim)
321     write(*,*)" ichim = ",ichim
322
323! Soil model
324     write(*,*)"Number of sub-surface layers for soil scheme?"
325     ! default value of nsoilmx set in comsoil_h
326     call getin_p("nsoilmx",nsoilmx)
327     write(*,*)" nsoilmx=",nsoilmx
328     
329     write(*,*)"Thickness of topmost soil layer (m)?"
330     ! default value of lay1_soil set in comsoil_h
331     call getin_p("lay1_soil",lay1_soil)
332     write(*,*)" lay1_soil = ",lay1_soil
333     
334     write(*,*)"Coefficient for soil layer thickness distribution?"
335     ! default value of alpha_soil set in comsoil_h
336     call getin_p("alpha_soil",alpha_soil)
337     write(*,*)" alpha_soil = ",alpha_soil
338
339     write(*,*)"Remove lower boundary?"
340     nosurf=.false.
341     call getin_p("nosurf",nosurf)
342     write(*,*)" nosurf = ",nosurf
343
344! Tests of incompatibility:
345     if (nosurf.and.callsoil) then
346       print*,'nosurf not compatible with soil scheme!'
347       print*,'... got to make a choice!'
348       call abort
349     endif
350
351     write(*,*)"Add an internal heat flux?", &
352                   "... matters only if callsoil=F"
353     intheat=0.
354     call getin_p("intheat",intheat)
355     write(*,*)" intheat = ",intheat
356
357     write(*,*)"Use Newtonian cooling for radiative transfer?"
358     newtonian=.false.
359     call getin_p("newtonian",newtonian)
360     write(*,*)" newtonian = ",newtonian
361
362! Tests of incompatibility:
363     if (newtonian.and.corrk) then
364       print*,'newtonian not compatible with correlated-k!'
365       call abort
366     endif
367     if (newtonian.and.calladj) then
368       print*,'newtonian not compatible with adjustment!'
369       call abort
370     endif
371     if (newtonian.and.calldifv) then
372       print*,'newtonian not compatible with a boundary layer!'
373       call abort
374     endif
375
376     write(*,*)"Test physics timescale in 1D?"
377     testradtimes=.false.
378     call getin_p("testradtimes",testradtimes)
379     write(*,*)" testradtimes = ",testradtimes
380
381! Test of incompatibility:
382! if testradtimes used, we must be in 1D
383     if (testradtimes.and.(ngrid.gt.1)) then
384       print*,'testradtimes can only be used in 1D!'
385       call abort
386     endif
387
388     write(*,*)"Default planetary temperature?"
389     tplanet=215.0
390     call getin_p("tplanet",tplanet)
391     write(*,*)" tplanet = ",tplanet
392
393     write(*,*)"Which star?"
394     startype=1 ! default value = Sol
395     call getin_p("startype",startype)
396     write(*,*)" startype = ",startype
397
398     write(*,*)"Value of stellar flux at 1 AU?"
399     Fat1AU=1356.0 ! default value = Sol today
400     call getin_p("Fat1AU",Fat1AU)
401     write(*,*)" Fat1AU = ",Fat1AU
402
403
404! TRACERS:
405
406!         write(*,*)"Number of radiatively active aerosols:"
407!         naerkind=0. ! default value
408!         call getin_p("naerkind",naerkind)
409!         write(*,*)" naerkind = ",naerkind
410
411 
412!=================================
413
414     write(*,*)"Radiatively active two-layer aersols?"
415     aeroback2lay=.false.     ! default value
416     call getin_p("aeroback2lay",aeroback2lay)
417     write(*,*)" aeroback2lay = ",aeroback2lay
418
419     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
420                    "in the tropospheric layer (visible)"
421     obs_tau_col_tropo=8.D0
422     call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
423     write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo
424
425     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
426                    "in the stratospheric layer (visible)"
427     obs_tau_col_strato=0.08D0
428     call getin_p("obs_tau_col_strato",obs_tau_col_strato)
429     write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato
430
431     write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
432     pres_bottom_tropo=66000.0
433     call getin_p("pres_bottom_tropo",pres_bottom_tropo)
434     write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo
435
436     write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa"
437     pres_top_tropo=18000.0
438     call getin_p("pres_top_tropo",pres_top_tropo)
439     write(*,*)" pres_top_tropo = ",pres_top_tropo
440
441     write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
442     pres_bottom_strato=2000.0
443     call getin_p("pres_bottom_strato",pres_bottom_strato)
444     write(*,*)" pres_bottom_strato = ",pres_bottom_strato
445
446     write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa"
447     pres_top_strato=100.0
448     call getin_p("pres_top_strato",pres_top_strato)
449     write(*,*)" pres_top_strato = ",pres_top_strato
450
451     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
452                    "tropospheric layer, in meters"
453     size_tropo=2.e-6
454     call getin_p("size_tropo",size_tropo)
455     write(*,*)" size_tropo = ",size_tropo
456
457     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
458                    "stratospheric layer, in meters"
459     size_strato=1.e-7
460     call getin_p("size_strato",size_strato)
461     write(*,*)" size_strato = ",size_strato
462
463!=================================
464
465     write(*,*) "Gravitationnal sedimentation ?"
466     sedimentation=.false. ! default value
467     call getin_p("sedimentation",sedimentation)
468     write(*,*) " sedimentation = ",sedimentation       
469
470     write(*,*) "Does user want to force cpp and mugaz?"
471     force_cpp=.false. ! default value
472     call getin_p("force_cpp",force_cpp)
473     write(*,*) " force_cpp = ",force_cpp
474
475     if (force_cpp) then
476       mugaz = -99999.
477       PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
478       call getin_p("mugaz",mugaz)
479       IF (mugaz.eq.-99999.) THEN
480           PRINT *, "mugaz must be set if force_cpp = T"
481           STOP
482       ELSE
483           write(*,*) "mugaz=",mugaz
484       ENDIF
485       cpp = -99999.
486       PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
487       call getin_p("cpp",cpp)
488       IF (cpp.eq.-99999.) THEN
489           PRINT *, "cpp must be set if force_cpp = T"
490           STOP
491       ELSE
492           write(*,*) "cpp=",cpp
493       ENDIF
494!         else
495!           mugaz=8.314*1000./pr
496     endif ! of if (force_cpp)
497     
498     
499     call su_gases(nlayer,tracer)     
500     call calc_cpp_mugaz
501     
502     
503     PRINT*,'--------------------------------------------'
504     PRINT*
505     PRINT*
506  ELSE
507     write(*,*)
508     write(*,*) 'Cannot read file callphys.def. Is it here ?'
509     stop
510  ENDIF ! of IF(iscallphys)
511
512  PRINT*
513  PRINT*,'inifis: daysec',daysec
514  PRINT*
515  PRINT*,'inifis: The radiative transfer is computed:'
516  PRINT*,'           each ',iradia,' physical time-step'
517  PRINT*,'        or each ',iradia*dtphys,' seconds'
518  PRINT*
519  PRINT*,'inifis: Chemistry is computed:'
520  PRINT*,'           each ',ichim,' physical time-step'
521  PRINT*,'        or each ',ichim*dtphys,' seconds'
522  PRINT*
523
524!-----------------------------------------------------------------------
525!     Some more initialization:
526!     ------------------------
527
528  ! Initializations for comgeomfi_h
529  totarea=SSUM(ngrid,parea,1)
530  call planetwide_sumval(parea,totarea_planet)
531
532  !! those are defined in comdiurn_h.F90
533  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
534  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
535  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
536  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
537
538  DO ig=1,ngrid
539     sinlat(ig)=sin(plat(ig))
540     coslat(ig)=cos(plat(ig))
541     sinlon(ig)=sin(plon(ig))
542     coslon(ig)=cos(plon(ig))
543  ENDDO
544
545  call ini_radinc_h(nlayer)
546 
547  ! allocate "radcommon_h" arrays
548  call ini_radcommon_h()
549
550  ! allocate "comsoil_h" arrays
551  call ini_comsoil_h(ngrid)
552     
553  END SUBROUTINE inifis
554
555END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.