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

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

+ Major clean of the new LMDZ.TITAN from too-generic options and routines (water, co2, ocean, surface type ...)
+ From this revision LMDZ.TITAN begins to be really separated from LMDZ.GENERIC
+ Partial desactivation of aerosols, only the dummy case is still enabled to keep the code running ( new aerosol routines to come in followings commits )

JVO

File size: 16.9 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
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
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! --------------------------------------------------------------
98!  Reading the "callphys.def" file controlling some key options
99! --------------------------------------------------------------
100     
101  ! check that 'callphys.def' file is around
102  OPEN(99,file='callphys.def',status='old',form='formatted',iostat=ierr)
103  CLOSE(99)
104  IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module
105     
106!!!      IF(ierr.EQ.0) THEN
107  IF(iscallphys) THEN
108     PRINT*
109     PRINT*
110     PRINT*,'--------------------------------------------'
111     PRINT*,' inifis: Parametres pour la physique (callphys.def)'
112     PRINT*,'--------------------------------------------'
113
114     write(*,*) "Directory where external input files are:"
115     ! default 'datadir' is set in "datadir_mod"
116     call getin_p("datadir",datadir) ! default path
117     write(*,*) " datadir = ",trim(datadir)
118
119     write(*,*) "Run with or without tracer transport ?"
120     tracer=.false. ! default value
121     call getin_p("tracer",tracer)
122     write(*,*) " tracer = ",tracer
123
124     write(*,*) "Run with or without atm mass update ", &
125            " due to tracer evaporation/condensation?"
126     mass_redistrib=.false. ! default value
127     call getin_p("mass_redistrib",mass_redistrib)
128     write(*,*) " mass_redistrib = ",mass_redistrib
129
130     write(*,*) "Diurnal cycle ?"
131     write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
132     diurnal=.true. ! default value
133     call getin_p("diurnal",diurnal)
134     write(*,*) " diurnal = ",diurnal
135
136     write(*,*) "Seasonal cycle ?"
137     write(*,*) "(if season=false, Ls stays constant, to value ", &
138         "set in 'start'"
139     season=.true. ! default value
140     call getin_p("season",season)
141     write(*,*) " season = ",season
142
143     write(*,*) "Tidally resonant rotation ?"
144     tlocked=.false. ! default value
145     call getin_p("tlocked",tlocked)
146     write(*,*) "tlocked = ",tlocked
147
148     write(*,*) "Saturn ring shadowing ?"
149     rings_shadow = .false.
150     call getin_p("rings_shadow", rings_shadow)
151     write(*,*) "rings_shadow = ", rings_shadow
152         
153     write(*,*) "Compute latitude-dependent gravity field?"
154     oblate = .false.
155     call getin_p("oblate", oblate)
156     write(*,*) "oblate = ", oblate
157
158     write(*,*) "Flattening of the planet (a-b)/a "
159     flatten = 0.0
160     call getin_p("flatten", flatten)
161     write(*,*) "flatten = ", flatten
162         
163
164     write(*,*) "Needed if oblate=.true.: J2"
165     J2 = 0.0
166     call getin_p("J2", J2)
167     write(*,*) "J2 = ", J2
168         
169     write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
170     MassPlanet = 0.0
171     call getin_p("MassPlanet", MassPlanet)
172     write(*,*) "MassPlanet = ", MassPlanet         
173
174     write(*,*) "Needed if oblate=.true.: Planet mean radius (m)"
175     Rmean = 0.0
176     call getin_p("Rmean", Rmean)
177     write(*,*) "Rmean = ", Rmean
178         
179! Test of incompatibility:
180! if tlocked, then diurnal should be false
181     if (tlocked.and.diurnal) then
182       print*,'If diurnal=true, we should turn off tlocked.'
183       stop
184     endif
185
186     write(*,*) "Tidal resonance ratio ?"
187     nres=0          ! default value
188     call getin_p("nres",nres)
189     write(*,*) "nres = ",nres
190
191     write(*,*) "Write some extra output to the screen ?"
192     lwrite=.false. ! default value
193     call getin_p("lwrite",lwrite)
194     write(*,*) " lwrite = ",lwrite
195
196     write(*,*) "Save statistics in file stats.nc ?"
197     callstats=.true. ! default value
198     call getin_p("callstats",callstats)
199     write(*,*) " callstats = ",callstats
200
201     write(*,*) "Test energy conservation of model physics ?"
202     enertest=.false. ! default value
203     call getin_p("enertest",enertest)
204     write(*,*) " enertest = ",enertest
205
206     write(*,*) "Check to see if cpp values used match gases.def ?"
207     check_cpp_match=.true. ! default value
208     call getin_p("check_cpp_match",check_cpp_match)
209     write(*,*) " check_cpp_match = ",check_cpp_match
210
211     write(*,*) "call radiative transfer ?"
212     callrad=.true. ! default value
213     call getin_p("callrad",callrad)
214     write(*,*) " callrad = ",callrad
215
216     write(*,*) "call correlated-k radiative transfer ?"
217     corrk=.true. ! default value
218     call getin_p("corrk",corrk)
219     write(*,*) " corrk = ",corrk
220
221     write(*,*) "prohibit calculations outside corrk T grid?"
222     strictboundcorrk=.true. ! default value
223     call getin_p("strictboundcorrk",strictboundcorrk)
224     write(*,*) "strictboundcorrk = ",strictboundcorrk
225
226     write(*,*) "call gaseous absorption in the visible bands?", &
227                    "(matters only if callrad=T)"
228     callgasvis=.false. ! default value
229     call getin_p("callgasvis",callgasvis)
230     write(*,*) " callgasvis = ",callgasvis
231       
232     write(*,*) "call continuum opacities in radiative transfer ?", &
233                    "(matters only if callrad=T)"
234     continuum=.true. ! default value
235     call getin_p("continuum",continuum)
236     write(*,*) " continuum = ",continuum
237 
238     write(*,*) "call turbulent vertical diffusion ?"
239     calldifv=.true. ! default value
240     call getin_p("calldifv",calldifv)
241     write(*,*) " calldifv = ",calldifv
242
243     write(*,*) "use turbdiff instead of vdifc ?"
244     UseTurbDiff=.true. ! default value
245     call getin_p("UseTurbDiff",UseTurbDiff)
246     write(*,*) " UseTurbDiff = ",UseTurbDiff
247
248     write(*,*) "call convective adjustment ?"
249     calladj=.true. ! default value
250     call getin_p("calladj",calladj)
251     write(*,*) " calladj = ",calladj
252
253     write(*,*) "Radiative timescale for Newtonian cooling ?"
254     tau_relax=30. ! default value
255     call getin_p("tau_relax",tau_relax)
256     write(*,*) " tau_relax = ",tau_relax
257     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
258
259     write(*,*)"call thermal conduction in the soil ?"
260     callsoil=.true. ! default value
261     call getin_p("callsoil",callsoil)
262     write(*,*) " callsoil = ",callsoil
263         
264     write(*,*)"Rad transfer is computed every iradia", &
265                   " physical timestep"
266     iradia=1 ! default value
267     call getin_p("iradia",iradia)
268     write(*,*)" iradia = ",iradia
269       
270     write(*,*)"Rayleigh scattering ?"
271     rayleigh=.false.
272     call getin_p("rayleigh",rayleigh)
273     write(*,*)" rayleigh = ",rayleigh
274
275     write(*,*) "Use blackbody for stellar spectrum ?"
276     stelbbody=.false. ! default value
277     call getin_p("stelbbody",stelbbody)
278     write(*,*) " stelbbody = ",stelbbody
279
280     write(*,*) "Stellar blackbody temperature ?"
281     stelTbb=5800.0 ! default value
282     call getin_p("stelTbb",stelTbb)
283     write(*,*) " stelTbb = ",stelTbb
284
285     write(*,*)"Output mean OLR in 1D?"
286     meanOLR=.false.
287     call getin_p("meanOLR",meanOLR)
288     write(*,*)" meanOLR = ",meanOLR
289
290     write(*,*)"Output spectral OLR in 3D?"
291     specOLR=.false.
292     call getin_p("specOLR",specOLR)
293     write(*,*)" specOLR = ",specOLR
294
295     write(*,*)"Uniform absorption in radiative transfer?"
296     graybody=.false.
297     call getin_p("graybody",graybody)
298     write(*,*)" graybody = ",graybody
299
300! Soil model
301     write(*,*)"Number of sub-surface layers for soil scheme?"
302     ! default value of nsoilmx set in comsoil_h
303     call getin_p("nsoilmx",nsoilmx)
304     write(*,*)" nsoilmx=",nsoilmx
305     
306     write(*,*)"Thickness of topmost soil layer (m)?"
307     ! default value of lay1_soil set in comsoil_h
308     call getin_p("lay1_soil",lay1_soil)
309     write(*,*)" lay1_soil = ",lay1_soil
310     
311     write(*,*)"Coefficient for soil layer thickness distribution?"
312     ! default value of alpha_soil set in comsoil_h
313     call getin_p("alpha_soil",alpha_soil)
314     write(*,*)" alpha_soil = ",alpha_soil
315
316     write(*,*)"Remove lower boundary?"
317     nosurf=.false.
318     call getin_p("nosurf",nosurf)
319     write(*,*)" nosurf = ",nosurf
320
321! Tests of incompatibility:
322     if (nosurf.and.callsoil) then
323       print*,'nosurf not compatible with soil scheme!'
324       print*,'... got to make a choice!'
325       call abort
326     endif
327
328     write(*,*)"Add an internal heat flux?", &
329                   "... matters only if callsoil=F"
330     intheat=0.
331     call getin_p("intheat",intheat)
332     write(*,*)" intheat = ",intheat
333
334     write(*,*)"Use Newtonian cooling for radiative transfer?"
335     newtonian=.false.
336     call getin_p("newtonian",newtonian)
337     write(*,*)" newtonian = ",newtonian
338
339! Tests of incompatibility:
340     if (newtonian.and.corrk) then
341       print*,'newtonian not compatible with correlated-k!'
342       call abort
343     endif
344     if (newtonian.and.calladj) then
345       print*,'newtonian not compatible with adjustment!'
346       call abort
347     endif
348     if (newtonian.and.calldifv) then
349       print*,'newtonian not compatible with a boundary layer!'
350       call abort
351     endif
352
353     write(*,*)"Test physics timescale in 1D?"
354     testradtimes=.false.
355     call getin_p("testradtimes",testradtimes)
356     write(*,*)" testradtimes = ",testradtimes
357
358! Test of incompatibility:
359! if testradtimes used, we must be in 1D
360     if (testradtimes.and.(ngrid.gt.1)) then
361       print*,'testradtimes can only be used in 1D!'
362       call abort
363     endif
364
365     write(*,*)"Default planetary temperature?"
366     tplanet=215.0
367     call getin_p("tplanet",tplanet)
368     write(*,*)" tplanet = ",tplanet
369
370     write(*,*)"Which star?"
371     startype=1 ! default value = Sol
372     call getin_p("startype",startype)
373     write(*,*)" startype = ",startype
374
375     write(*,*)"Value of stellar flux at 1 AU?"
376     Fat1AU=1356.0 ! default value = Sol today
377     call getin_p("Fat1AU",Fat1AU)
378     write(*,*)" Fat1AU = ",Fat1AU
379
380
381! TRACERS:
382
383!         write(*,*)"Number of radiatively active aerosols:"
384!         naerkind=0. ! default value
385!         call getin_p("naerkind",naerkind)
386!         write(*,*)" naerkind = ",naerkind
387
388 
389!=================================
390
391     write(*,*)"Radiatively active two-layer aersols?"
392     aeroback2lay=.false.     ! default value
393     call getin_p("aeroback2lay",aeroback2lay)
394     write(*,*)" aeroback2lay = ",aeroback2lay
395
396     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
397                    "in the tropospheric layer (visible)"
398     obs_tau_col_tropo=8.D0
399     call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
400     write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo
401
402     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
403                    "in the stratospheric layer (visible)"
404     obs_tau_col_strato=0.08D0
405     call getin_p("obs_tau_col_strato",obs_tau_col_strato)
406     write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato
407
408     write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
409     pres_bottom_tropo=66000.0
410     call getin_p("pres_bottom_tropo",pres_bottom_tropo)
411     write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo
412
413     write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa"
414     pres_top_tropo=18000.0
415     call getin_p("pres_top_tropo",pres_top_tropo)
416     write(*,*)" pres_top_tropo = ",pres_top_tropo
417
418     write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
419     pres_bottom_strato=2000.0
420     call getin_p("pres_bottom_strato",pres_bottom_strato)
421     write(*,*)" pres_bottom_strato = ",pres_bottom_strato
422
423     write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa"
424     pres_top_strato=100.0
425     call getin_p("pres_top_strato",pres_top_strato)
426     write(*,*)" pres_top_strato = ",pres_top_strato
427
428     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
429                    "tropospheric layer, in meters"
430     size_tropo=2.e-6
431     call getin_p("size_tropo",size_tropo)
432     write(*,*)" size_tropo = ",size_tropo
433
434     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
435                    "stratospheric layer, in meters"
436     size_strato=1.e-7
437     call getin_p("size_strato",size_strato)
438     write(*,*)" size_strato = ",size_strato
439
440!=================================
441
442     write(*,*)"What is the saturation % of the variable species?"
443     satval=0.8
444     call getin_p("satval",satval)
445     write(*,*)" satval = ",satval
446
447     write(*,*) "Gravitationnal sedimentation ?"
448     sedimentation=.false. ! default value
449     call getin_p("sedimentation",sedimentation)
450     write(*,*) " sedimentation = ",sedimentation       
451
452     write(*,*) "Does user want to force cpp and mugaz?"
453     force_cpp=.false. ! default value
454     call getin_p("force_cpp",force_cpp)
455     write(*,*) " force_cpp = ",force_cpp
456
457     if (force_cpp) then
458       mugaz = -99999.
459       PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
460       call getin_p("mugaz",mugaz)
461       IF (mugaz.eq.-99999.) THEN
462           PRINT *, "mugaz must be set if force_cpp = T"
463           STOP
464       ELSE
465           write(*,*) "mugaz=",mugaz
466       ENDIF
467       cpp = -99999.
468       PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
469       call getin_p("cpp",cpp)
470       IF (cpp.eq.-99999.) THEN
471           PRINT *, "cpp must be set if force_cpp = T"
472           STOP
473       ELSE
474           write(*,*) "cpp=",cpp
475       ENDIF
476!         else
477!           mugaz=8.314*1000./pr
478     endif ! of if (force_cpp)
479     call su_gases
480     call calc_cpp_mugaz
481
482     PRINT*,'--------------------------------------------'
483     PRINT*
484     PRINT*
485  ELSE
486     write(*,*)
487     write(*,*) 'Cannot read file callphys.def. Is it here ?'
488     stop
489  ENDIF ! of IF(iscallphys)
490
491  PRINT*
492  PRINT*,'inifis: daysec',daysec
493  PRINT*
494  PRINT*,'inifis: The radiative transfer is computed:'
495  PRINT*,'           each ',iradia,' physical time-step'
496  PRINT*,'        or each ',iradia*dtphys,' seconds'
497  PRINT*
498
499
500!-----------------------------------------------------------------------
501!     Some more initialization:
502!     ------------------------
503
504  ! Initializations for comgeomfi_h
505  totarea=SSUM(ngrid,parea,1)
506  call planetwide_sumval(parea,totarea_planet)
507
508  !! those are defined in comdiurn_h.F90
509  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
510  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
511  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
512  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
513
514  DO ig=1,ngrid
515     sinlat(ig)=sin(plat(ig))
516     coslat(ig)=cos(plat(ig))
517     sinlon(ig)=sin(plon(ig))
518     coslon(ig)=cos(plon(ig))
519  ENDDO
520
521  ! initialize variables in radinc_h
522  call ini_radinc_h(nlayer)
523 
524  ! allocate "radcommon_h" arrays
525  call ini_radcommon_h()
526
527  ! allocate "comsoil_h" arrays
528  call ini_comsoil_h(ngrid)
529     
530  END SUBROUTINE inifis
531
532END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.