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

Last change on this file since 1788 was 1788, checked in by jvatant, 7 years ago

Get rid of all the old-generic dummy aerosol scheme ( just left scatterers_h for compilation )
as the new microphysics for Titan will be plugged in
-> even removed sedimentation ( will be done in the microphysical model )
--JVO

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