1 | MODULE inifis_mod |
---|
2 | IMPLICIT NONE |
---|
3 | |
---|
4 | CONTAINS |
---|
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 | |
---|
480 | END MODULE inifis_mod |
---|