1 | MODULE conf_phys_mod |
---|
2 | |
---|
3 | IMPLICIT NONE |
---|
4 | |
---|
5 | CONTAINS |
---|
6 | |
---|
7 | SUBROUTINE conf_phys(ngrid,nlayer,nq) |
---|
8 | |
---|
9 | !======================================================================= |
---|
10 | ! |
---|
11 | ! purpose: |
---|
12 | ! ------- |
---|
13 | ! |
---|
14 | ! Initialisation for the physical parametrisations |
---|
15 | ! flags (i.e. run-time options) of the Mars PCM |
---|
16 | !----------------------------------------------------------------------- |
---|
17 | |
---|
18 | USE ioipsl_getin_p_mod, ONLY : getin_p |
---|
19 | use tracer_mod, only : nuice_sed, ccn_factor, nuiceco2_sed, |
---|
20 | & nuice_ref,nuiceco2_ref |
---|
21 | use surfdat_h, only: albedo_h2o_cap,albedo_h2o_frost, |
---|
22 | & frost_albedo_threshold, inert_h2o_ice, |
---|
23 | & frost_metam_threshold,old_wsublimation_scheme |
---|
24 | use time_phylmdz_mod, only: ecritphy,day_step,iphysiq,ecritstart, |
---|
25 | & daysec,dtphys |
---|
26 | use dimradmars_mod, only: naerkind, name_iaer, |
---|
27 | & ini_scatterers,tauvis |
---|
28 | use datafile_mod, only: datadir |
---|
29 | use wstats_mod, only: callstats |
---|
30 | use writediagsoil_mod, only: diagsoil |
---|
31 | use calchim_mod, only: ichemistry |
---|
32 | use co2condens_mod, only: scavco2cond |
---|
33 | use dust_param_mod, only: dustbin, doubleq, submicron, active, |
---|
34 | & lifting, freedust, callddevil, |
---|
35 | & dustscaling_mode, |
---|
36 | & reff_driven_IRtoVIS_scenario |
---|
37 | use aeropacity_mod, only: iddist, topdustref |
---|
38 | USE mod_phys_lmdz_transfert_para, ONLY: bcast |
---|
39 | USE paleoclimate_mod,ONLY: paleoclimate,albedo_perenialco2, |
---|
40 | & lag_layer |
---|
41 | use microphys_h, only: mteta |
---|
42 | use comsoil_h, only: adsorption_soil, choice_ads |
---|
43 | |
---|
44 | IMPLICIT NONE |
---|
45 | |
---|
46 | include "callkeys.h" |
---|
47 | |
---|
48 | INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns |
---|
49 | INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers |
---|
50 | INTEGER,INTENT(IN) :: nq ! number of tracers |
---|
51 | |
---|
52 | INTEGER ierr,j |
---|
53 | character(len=20),parameter :: modname="conf_phys" |
---|
54 | |
---|
55 | CHARACTER ch1*12 |
---|
56 | #ifndef MESOSCALE |
---|
57 | ! read in some parameters from "run.def" for physics, |
---|
58 | ! or shared between dynamics and physics. |
---|
59 | ecritphy=240 ! default value |
---|
60 | call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics, |
---|
61 | ! in dynamical steps |
---|
62 | day_step=960 ! default value |
---|
63 | call getin_p("day_step",day_step) ! number of dynamical steps per day |
---|
64 | iphysiq=20 ! default value |
---|
65 | call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step |
---|
66 | ecritstart=0 ! default value |
---|
67 | call getin_p("ecritstart",ecritstart) ! write a restart every ecristart steps |
---|
68 | #endif |
---|
69 | |
---|
70 | ! -------------------------------------------------------------- |
---|
71 | ! Reading the "callphys.def" file controlling some key options |
---|
72 | ! -------------------------------------------------------------- |
---|
73 | |
---|
74 | !$OMP MASTER |
---|
75 | ! check that 'callphys.def' file is around |
---|
76 | OPEN(99,file='callphys.def',status='old',form='formatted' |
---|
77 | & ,iostat=ierr) |
---|
78 | CLOSE(99) |
---|
79 | !$OMP END MASTER |
---|
80 | call bcast(ierr) |
---|
81 | ! ierr=0 |
---|
82 | |
---|
83 | IF(ierr.EQ.0) THEN |
---|
84 | PRINT* |
---|
85 | PRINT* |
---|
86 | PRINT*,'--------------------------------------------' |
---|
87 | PRINT*,' conf_phys: Parameters for the physics (callphys.def)' |
---|
88 | PRINT*,'--------------------------------------------' |
---|
89 | |
---|
90 | write(*,*) "Directory where external input files are:" |
---|
91 | ! default path is set in datafile_mod |
---|
92 | call getin_p("datadir",datadir) |
---|
93 | write(*,*) " datadir = ",trim(datadir) |
---|
94 | |
---|
95 | write(*,*) "Initialize physics with startfi.nc file ?" |
---|
96 | startphy_file=.true. |
---|
97 | call getin_p("startphy_file",startphy_file) |
---|
98 | write(*,*) "startphy_file", startphy_file |
---|
99 | |
---|
100 | write(*,*) "Diurnal cycle ?" |
---|
101 | write(*,*) "(if diurnal=False, diurnal averaged solar heating)" |
---|
102 | diurnal=.true. ! default value |
---|
103 | call getin_p("diurnal",diurnal) |
---|
104 | write(*,*) " diurnal = ",diurnal |
---|
105 | |
---|
106 | write(*,*) "Seasonal cycle ?" |
---|
107 | write(*,*) "(if season=False, Ls stays constant, to value ", |
---|
108 | & "set in 'start'" |
---|
109 | season=.true. ! default value |
---|
110 | call getin_p("season",season) |
---|
111 | write(*,*) " season = ",season |
---|
112 | |
---|
113 | write(*,*) "Write some extra output to the screen ?" |
---|
114 | lwrite=.false. ! default value |
---|
115 | call getin_p("lwrite",lwrite) |
---|
116 | write(*,*) " lwrite = ",lwrite |
---|
117 | |
---|
118 | write(*,*) "Save statistics in file stats.nc ?" |
---|
119 | #ifdef MESOSCALE |
---|
120 | callstats=.false. ! default value |
---|
121 | #else |
---|
122 | callstats=.true. ! default value |
---|
123 | #endif |
---|
124 | call getin_p("callstats",callstats) |
---|
125 | write(*,*) " callstats = ",callstats |
---|
126 | |
---|
127 | write(*,*) "Write sub-surface fields in file diagsoil.nc ?" |
---|
128 | diagsoil=.false. ! default value |
---|
129 | call getin_p("diagsoil",diagsoil) |
---|
130 | write(*,*) " diagsoil = ",diagsoil |
---|
131 | |
---|
132 | write(*,*) "Save EOF profiles in file 'profiles' for ", |
---|
133 | & "Climate Database?" |
---|
134 | calleofdump=.false. ! default value |
---|
135 | call getin_p("calleofdump",calleofdump) |
---|
136 | write(*,*) " calleofdump = ",calleofdump |
---|
137 | |
---|
138 | write(*,*) "Dust scenario: 1=constant dust (read from startfi", |
---|
139 | & " or set as tauvis); 2=Viking scenario; =3 MGS scenario,", |
---|
140 | & "=6 cold (low dust) scenario; =7 warm (high dust) scenario ", |
---|
141 | & "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation" |
---|
142 | iaervar=3 ! default value |
---|
143 | call getin_p("iaervar",iaervar) |
---|
144 | write(*,*) " iaervar = ",iaervar |
---|
145 | |
---|
146 | write(*,*) "Reference (visible) dust opacity at 610 Pa ", |
---|
147 | & "(matters only if iaervar=1)" |
---|
148 | ! NB: default value of tauvis is set/read in startfi.nc file |
---|
149 | call getin_p("tauvis",tauvis) |
---|
150 | write(*,*) " tauvis = ",tauvis |
---|
151 | |
---|
152 | write(*,*) "Dust vertical distribution:" |
---|
153 | write(*,*) "(=1 top set by topdustref parameter;", |
---|
154 | & " =2 Viking scenario; =3 MGS scenario)" |
---|
155 | iddist=3 ! default value |
---|
156 | call getin_p("iddist",iddist) |
---|
157 | write(*,*) " iddist = ",iddist |
---|
158 | |
---|
159 | write(*,*) "Dust top altitude (km). (Matters only if iddist=1)" |
---|
160 | topdustref= 90.0 ! default value |
---|
161 | call getin_p("topdustref",topdustref) |
---|
162 | write(*,*) " topdustref = ",topdustref |
---|
163 | |
---|
164 | write(*,*) "Prescribed surface thermal flux (H/(rho*cp),K m/s)" |
---|
165 | tke_heat_flux=0. ! default value |
---|
166 | call getin_p("tke_heat_flux",tke_heat_flux) |
---|
167 | write(*,*) " tke_heat_flux = ",tke_heat_flux |
---|
168 | write(*,*) " 0 means the usual schemes are computing" |
---|
169 | |
---|
170 | write(*,*) "call radiative transfer ?" |
---|
171 | callrad=.true. ! default value |
---|
172 | call getin_p("callrad",callrad) |
---|
173 | write(*,*) " callrad = ",callrad |
---|
174 | |
---|
175 | write(*,*) "call slope insolation scheme ?", |
---|
176 | & "(matters only if callrad=T)" |
---|
177 | #ifdef MESOSCALE |
---|
178 | callslope=.true. ! default value |
---|
179 | #else |
---|
180 | callslope=.false. ! default value (not supported yet) |
---|
181 | #endif |
---|
182 | call getin_p("callslope",callslope) |
---|
183 | write(*,*) " callslope = ",callslope |
---|
184 | |
---|
185 | write(*,*) "call NLTE radiative schemes ?", |
---|
186 | & "(matters only if callrad=T)" |
---|
187 | callnlte=.false. ! default value |
---|
188 | call getin_p("callnlte",callnlte) |
---|
189 | write(*,*) " callnlte = ",callnlte |
---|
190 | |
---|
191 | nltemodel=0 !default value |
---|
192 | write(*,*) "NLTE model?" |
---|
193 | write(*,*) "0 -> old model, static O" |
---|
194 | write(*,*) "1 -> old model, dynamic O" |
---|
195 | write(*,*) "2 -> new model" |
---|
196 | write(*,*) "(matters only if callnlte=T)" |
---|
197 | call getin_p("nltemodel",nltemodel) |
---|
198 | write(*,*) " nltemodel = ",nltemodel |
---|
199 | |
---|
200 | write(*,*) "call CO2 NIR absorption ?", |
---|
201 | & "(matters only if callrad=T)" |
---|
202 | callnirco2=.false. ! default value |
---|
203 | call getin_p("callnirco2",callnirco2) |
---|
204 | write(*,*) " callnirco2 = ",callnirco2 |
---|
205 | |
---|
206 | write(*,*) "New NIR NLTE correction ?", |
---|
207 | $ "0-> old model (no correction)", |
---|
208 | $ "1-> new correction", |
---|
209 | $ "(matters only if callnirco2=T)" |
---|
210 | #ifdef MESOSCALE |
---|
211 | nircorr=0 !default value. this is OK below 60 km. |
---|
212 | #else |
---|
213 | nircorr=0 !default value |
---|
214 | #endif |
---|
215 | call getin_p("nircorr",nircorr) |
---|
216 | write(*,*) " nircorr = ",nircorr |
---|
217 | |
---|
218 | write(*,*) "call turbulent vertical diffusion ?" |
---|
219 | calldifv=.true. ! default value |
---|
220 | call getin_p("calldifv",calldifv) |
---|
221 | write(*,*) " calldifv = ",calldifv |
---|
222 | |
---|
223 | write(*,*) "call thermals ?" |
---|
224 | calltherm=.false. ! default value |
---|
225 | call getin_p("calltherm",calltherm) |
---|
226 | write(*,*) " calltherm = ",calltherm |
---|
227 | |
---|
228 | write(*,*) "call convective adjustment ?" |
---|
229 | calladj=.true. ! default value |
---|
230 | call getin_p("calladj",calladj) |
---|
231 | write(*,*) " calladj = ",calladj |
---|
232 | |
---|
233 | if (calltherm .and. calladj) then |
---|
234 | print*,'!!! PLEASE NOTE !!!' |
---|
235 | print*,'convective adjustment is on' |
---|
236 | print*,'but since thermal plume model is on' |
---|
237 | print*,'convadj is only activated above the PBL' |
---|
238 | endif |
---|
239 | |
---|
240 | write(*,*) "used latest version of yamada scheme?" |
---|
241 | callyamada4=.true. ! default value |
---|
242 | call getin_p("callyamada4",callyamada4) |
---|
243 | write(*,*) " callyamada4 = ",callyamada4 |
---|
244 | |
---|
245 | if (calltherm .and. .not.callyamada4) then |
---|
246 | print*,'!!!! WARNING WARNING WARNING !!!!' |
---|
247 | print*,'if calltherm=T we strongly advise that ' |
---|
248 | print*,'you set the flag callyamada4 to T ' |
---|
249 | print*,'!!!! WARNING WARNING WARNING !!!!' |
---|
250 | endif |
---|
251 | |
---|
252 | write(*,*) "call Richardson-based surface layer ?" |
---|
253 | callrichsl=.false. ! default value |
---|
254 | call getin_p("callrichsl",callrichsl) |
---|
255 | write(*,*) " callrichsl = ",callrichsl |
---|
256 | |
---|
257 | if (calltherm .and. .not.callrichsl) then |
---|
258 | print*,'WARNING WARNING WARNING' |
---|
259 | print*,'if calltherm=T we strongly advise that ' |
---|
260 | print*,'you use the new surface layer scheme ' |
---|
261 | print*,'by setting callrichsl=T ' |
---|
262 | endif |
---|
263 | |
---|
264 | if (calladj .and. callrichsl .and. (.not. calltherm)) then |
---|
265 | print*,'You should not be calling the convective adjustment |
---|
266 | & scheme with the Richardson surface-layer and without the thermals |
---|
267 | &. This approach is not |
---|
268 | & physically consistent and can lead to unrealistic friction |
---|
269 | & values.' |
---|
270 | print*,'If you want to use the Ri. surface-layer, either |
---|
271 | & activate thermals OR de-activate the convective adjustment.' |
---|
272 | call abort_physic(modname, |
---|
273 | & "Richardson layer must be used with thermals",1) |
---|
274 | endif |
---|
275 | |
---|
276 | write(*,*) "call CO2 condensation ?" |
---|
277 | callcond=.true. ! default value |
---|
278 | call getin_p("callcond",callcond) |
---|
279 | write(*,*) " callcond = ",callcond |
---|
280 | |
---|
281 | write(*,*)"call thermal conduction in the soil ?" |
---|
282 | callsoil=.true. ! default value |
---|
283 | call getin_p("callsoil",callsoil) |
---|
284 | write(*,*) " callsoil = ",callsoil |
---|
285 | |
---|
286 | |
---|
287 | write(*,*)"call Lott's gravity wave/subgrid topography ", |
---|
288 | & "scheme ?" |
---|
289 | calllott=.true. ! default value |
---|
290 | call getin_p("calllott",calllott) |
---|
291 | write(*,*)" calllott = ",calllott |
---|
292 | |
---|
293 | write(*,*)"call Lott's non-oro GWs parameterisation ", |
---|
294 | & "scheme ?" |
---|
295 | calllott_nonoro=.false. ! default value |
---|
296 | call getin_p("calllott_nonoro",calllott_nonoro) |
---|
297 | write(*,*)" calllott_nonoro = ",calllott_nonoro |
---|
298 | |
---|
299 | ! rocket dust storm injection scheme |
---|
300 | write(*,*)"call rocket dust storm parametrization" |
---|
301 | rdstorm=.false. ! default value |
---|
302 | call getin_p("rdstorm",rdstorm) |
---|
303 | write(*,*)" rdstorm = ",rdstorm |
---|
304 | ! rocket dust storm detrainment coefficient |
---|
305 | coeff_detrainment=0.02 ! default value |
---|
306 | call getin_p("coeff_detrainment",coeff_detrainment) |
---|
307 | write(*,*)" coeff_detrainment = ",coeff_detrainment |
---|
308 | |
---|
309 | ! entrainment by mountain top dust flows scheme |
---|
310 | write(*,*)"call mountain top dust flows parametrization" |
---|
311 | topflows=.false. ! default value |
---|
312 | call getin_p("topflows",topflows) |
---|
313 | write(*,*)" topflows = ",topflows |
---|
314 | |
---|
315 | ! latent heat release from ground water ice sublimation/condensation |
---|
316 | write(*,*)"latent heat release during sublimation", |
---|
317 | & " /condensation of ground water ice" |
---|
318 | latentheat_surfwater=.true. ! default value |
---|
319 | call getin_p("latentheat_surfwater",latentheat_surfwater) |
---|
320 | write(*,*)" latentheat_surfwater = ",latentheat_surfwater |
---|
321 | |
---|
322 | write(*,*)"rad.transfer is computed every iradia", |
---|
323 | & " physical timestep" |
---|
324 | iradia=1 ! default value |
---|
325 | call getin_p("iradia",iradia) |
---|
326 | write(*,*)" iradia = ",iradia |
---|
327 | |
---|
328 | |
---|
329 | write(*,*)"Output of the exchange coefficient mattrix ?", |
---|
330 | & "(for diagnostics only)" |
---|
331 | callg2d=.false. ! default value |
---|
332 | call getin_p("callg2d",callg2d) |
---|
333 | write(*,*)" callg2d = ",callg2d |
---|
334 | |
---|
335 | write(*,*)"Rayleigh scattering : (should be .false. for now)" |
---|
336 | rayleigh=.false. |
---|
337 | call getin_p("rayleigh",rayleigh) |
---|
338 | write(*,*)" rayleigh = ",rayleigh |
---|
339 | |
---|
340 | ! PALEOCLIMATE |
---|
341 | |
---|
342 | write(*,*)"Using lag layer??" |
---|
343 | lag_layer=.false. |
---|
344 | call getin_p("lag_layer",lag_layer) |
---|
345 | write(*,*) " lag_layer = ", lag_layer |
---|
346 | |
---|
347 | write(*,*)"Is it paleoclimate run?" |
---|
348 | paleoclimate=.false. ! default value |
---|
349 | call getin_p("paleoclimate",paleoclimate) |
---|
350 | write(*,*)" paleoclimate = ",paleoclimate |
---|
351 | |
---|
352 | |
---|
353 | write(*,*)"Albedo for perenial CO2 ice?" |
---|
354 | albedo_perenialco2 = 0.85 ! default value |
---|
355 | call getin_p("albedo_perenialco2",albedo_perenialco2) |
---|
356 | write(*,*)"albedo_perenialco2 = ",albedo_perenialco2 |
---|
357 | |
---|
358 | ! TRACERS: |
---|
359 | |
---|
360 | ! dustbin |
---|
361 | write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)" |
---|
362 | dustbin=0 ! default value |
---|
363 | call getin_p("dustbin",dustbin) |
---|
364 | write(*,*)" dustbin = ",dustbin |
---|
365 | ! active |
---|
366 | write(*,*)"Radiatively active dust ? (matters if dustbin>0)" |
---|
367 | active=.false. ! default value |
---|
368 | call getin_p("active",active) |
---|
369 | write(*,*)" active = ",active |
---|
370 | |
---|
371 | ! Test of incompatibility: |
---|
372 | ! if active is used, then dustbin should be > 0 |
---|
373 | |
---|
374 | if (active.and.(dustbin.lt.1)) then |
---|
375 | print*,'if active is used, then dustbin should > 0' |
---|
376 | call abort_physic(modname, |
---|
377 | & "active option requires dustbin < 0",1) |
---|
378 | endif |
---|
379 | ! doubleq |
---|
380 | write(*,*)"use mass and number mixing ratios to predict", |
---|
381 | & " dust size ?" |
---|
382 | doubleq=.false. ! default value |
---|
383 | call getin_p("doubleq",doubleq) |
---|
384 | write(*,*)" doubleq = ",doubleq |
---|
385 | ! submicron |
---|
386 | submicron=.false. ! default value |
---|
387 | call getin_p("submicron",submicron) |
---|
388 | write(*,*)" submicron = ",submicron |
---|
389 | |
---|
390 | ! Test of incompatibility: |
---|
391 | ! if doubleq is used, then dustbin should be 2 |
---|
392 | |
---|
393 | if (doubleq.and.(dustbin.ne.2)) then |
---|
394 | print*,'if doubleq is used, then dustbin should be 2' |
---|
395 | call abort_physic(modname, |
---|
396 | & "doubleq option requires dustbin = 2",1) |
---|
397 | endif |
---|
398 | if (doubleq.and.submicron.and.(nq.LT.3)) then |
---|
399 | print*,'If doubleq is used with a submicron tracer,' |
---|
400 | print*,' then the number of tracers has to be' |
---|
401 | print*,' larger than 3.' |
---|
402 | call abort_physic(modname, |
---|
403 | & "submicron option requires dustbin > 2",1) |
---|
404 | endif |
---|
405 | |
---|
406 | ! lifting |
---|
407 | write(*,*)"dust lifted by GCM surface winds ?" |
---|
408 | lifting=.false. ! default value |
---|
409 | call getin_p("lifting",lifting) |
---|
410 | write(*,*)" lifting = ",lifting |
---|
411 | |
---|
412 | ! Test of incompatibility: |
---|
413 | ! if lifting is used, then dustbin should be > 0 |
---|
414 | |
---|
415 | if (lifting.and.(dustbin.lt.1)) then |
---|
416 | print*,'if lifting is used, then dustbin should > 0' |
---|
417 | call abort_physic(modname, |
---|
418 | & "lifting option requires dustbin > 0",1) |
---|
419 | endif |
---|
420 | |
---|
421 | ! dust injection scheme |
---|
422 | dustinjection=0 ! default: no injection scheme |
---|
423 | call getin_p("dustinjection",dustinjection) |
---|
424 | write(*,*)" dustinjection = ",dustinjection |
---|
425 | ! dust injection scheme coefficient |
---|
426 | coeff_injection=0.25 ! default value |
---|
427 | call getin_p("coeff_injection",coeff_injection) |
---|
428 | write(*,*)" coeff_in,jection = ",coeff_injection |
---|
429 | ! timing for dust injection |
---|
430 | ti_injection=0. ! default value |
---|
431 | tf_injection=24. ! default value |
---|
432 | call getin_p("ti_injection",ti_injection) |
---|
433 | write(*,*)" ti_injection = ",ti_injection |
---|
434 | call getin_p("tf_injection",tf_injection) |
---|
435 | write(*,*)" tf_injection = ",tf_injection |
---|
436 | |
---|
437 | ! free evolving dust |
---|
438 | ! freedust=true just says that there is no lifting and no dust opacity scaling. |
---|
439 | write(*,*)"dust lifted by GCM surface winds ?" |
---|
440 | freedust=.false. ! default value |
---|
441 | call getin_p("freedust",freedust) |
---|
442 | write(*,*)" freedust = ",freedust |
---|
443 | if (freedust.and..not.doubleq) then |
---|
444 | print*,'freedust should be used with doubleq !' |
---|
445 | call abort_physic(modname, |
---|
446 | & "freedust option requires doubleq",1) |
---|
447 | endif |
---|
448 | |
---|
449 | ! dust rescaling mode (if any) |
---|
450 | ! =0 --> freedust, tauscaling=1 |
---|
451 | ! =1 --> GCM5.3-like, tauscaling computed wrt tau_pref_scenario |
---|
452 | ! =2 --> tauscaling=1 but dust_rad_adjust (rescaling for the radiative transfer only) |
---|
453 | if (freedust) then |
---|
454 | dustscaling_mode=0 |
---|
455 | else |
---|
456 | dustscaling_mode=1 ! GCMv5.3 style |
---|
457 | endif |
---|
458 | call getin_p("dustscaling_mode",dustscaling_mode) |
---|
459 | write(*,*) "dustscaling_mode=",dustscaling_mode |
---|
460 | if ((dustscaling_mode.eq.1).and.freedust) then |
---|
461 | print*,'freedust and dustscaling_mode=1 are incompatible !' |
---|
462 | call abort_physic(modname, |
---|
463 | & "freedust and dustscaling_mode=1 are incompatible",1) |
---|
464 | endif |
---|
465 | |
---|
466 | #ifndef MESOSCALE |
---|
467 | ! this test is valid in GCM case |
---|
468 | ! ... not in mesoscale case, for we want to activate mesoscale lifting |
---|
469 | if (freedust.and.dustinjection.eq.0)then |
---|
470 | if(lifting) then |
---|
471 | print*,'if freedust is used and dustinjection = 0, |
---|
472 | & then lifting should not be used' |
---|
473 | call abort_physic(modname, |
---|
474 | & "freedust option with dustinjection = 0"// |
---|
475 | & " requires lifting to be false",1) |
---|
476 | endif |
---|
477 | endif |
---|
478 | #endif |
---|
479 | if (dustinjection.eq.1)then |
---|
480 | if(.not.lifting) then |
---|
481 | print*,"if dustinjection=1, then lifting should be true" |
---|
482 | call abort_physic(modname, |
---|
483 | & "dustinjection=1 requires lifting",1) |
---|
484 | endif |
---|
485 | if(.not.freedust) then |
---|
486 | print*,"if dustinjection=1, then freedust should be true" |
---|
487 | call abort_physic(modname, |
---|
488 | & "dustinjection=1 requires freedust",1) |
---|
489 | endif |
---|
490 | endif |
---|
491 | ! rocket dust storm and entrainment by top flows |
---|
492 | ! Test of incompatibility: |
---|
493 | ! if rdstorm or topflows is used, then doubleq should be true |
---|
494 | if ((rdstorm.or.topflows).and..not.doubleq) then |
---|
495 | print*,'if rdstorm or topflows is used, then doubleq |
---|
496 | & should be used !' |
---|
497 | call abort_physic(modname, |
---|
498 | & "rdstorm or topflows requires doubleq",1) |
---|
499 | endif |
---|
500 | if ((rdstorm.or.topflows).and..not.active) then |
---|
501 | print*,'if rdstorm or topflows is used, then active |
---|
502 | & should be used !' |
---|
503 | call abort_physic(modname, |
---|
504 | & "rdstorm or topflows requires active",1) |
---|
505 | endif |
---|
506 | if (rdstorm.and..not.lifting) then |
---|
507 | print*,'if rdstorm is used, then lifting |
---|
508 | & should be used !' |
---|
509 | call abort_physic(modname, |
---|
510 | & "rdstorm requires lifting",1) |
---|
511 | endif |
---|
512 | if ((rdstorm.or.topflows).and..not.freedust) then |
---|
513 | print*,'if rdstorm or topflows is used, then freedust |
---|
514 | & should be used !' |
---|
515 | call abort_physic(modname, |
---|
516 | & "rdstorm or topflows requires freedust",1) |
---|
517 | endif |
---|
518 | if (rdstorm.and.(dustinjection.eq.0)) then |
---|
519 | print*,'if rdstorm is used, then dustinjection |
---|
520 | & should be used !' |
---|
521 | call abort_physic(modname, |
---|
522 | & "rdstorm requires dustinjection",1) |
---|
523 | endif |
---|
524 | ! Dust IR opacity |
---|
525 | write(*,*)" Wavelength for infrared opacity of dust ?" |
---|
526 | write(*,*)" Choices are:" |
---|
527 | write(*,*)" tes --- > 9.3 microns [default]" |
---|
528 | write(*,*)" mcs --- > 21.6 microns" |
---|
529 | ! |
---|
530 | ! WARNING WARNING WARNING WARNING WARNING WARNING |
---|
531 | ! |
---|
532 | ! BEFORE ADDING A NEW VALUE, BE SURE THAT THE |
---|
533 | ! CORRESPONDING WAVELENGTH IS IN THE LOOKUP TABLE, |
---|
534 | ! OR AT LEAST NO TO FAR, TO AVOID FALLACIOUS INTERPOLATIONS. |
---|
535 | ! |
---|
536 | dustiropacity="tes" !default value |
---|
537 | call getin_p("dustiropacity",dustiropacity) |
---|
538 | write(*,*)" dustiropacity = ",trim(dustiropacity) |
---|
539 | select case (trim(dustiropacity)) |
---|
540 | case ("tes") |
---|
541 | dustrefir = 9.3E-6 |
---|
542 | case ("mcs") |
---|
543 | dustrefir = 21.6E-6 |
---|
544 | case default |
---|
545 | write(*,*) trim(dustiropacity), |
---|
546 | & " is not a valid option for dustiropacity" |
---|
547 | call abort_physic(modname, |
---|
548 | & "invalid dustiropacity option value",1) |
---|
549 | end select |
---|
550 | ! Dust scenario IR to VIS conversion |
---|
551 | write(*,*)"Use an IR to VIS conversion coefficient" |
---|
552 | write(*,*)"for the dust scenario, that is dependent" |
---|
553 | write(*,*)"on the GCM dust effective radius," |
---|
554 | write(*,*)"instead of a fixed 2.6 coefficient ?" |
---|
555 | reff_driven_IRtoVIS_scenario=.false. !default value |
---|
556 | call getin_p("reff_driven_IRtoVIS_scenario", |
---|
557 | & reff_driven_IRtoVIS_scenario) |
---|
558 | write(*,*)" reff_driven_IRtoVIS_scenario = ", |
---|
559 | & reff_driven_IRtoVIS_scenario |
---|
560 | ! Test of incompatibility: |
---|
561 | ! if reff_driven_IRtoVIS_scenario=.true., |
---|
562 | ! dustrefir must be 9.3E-6 = scenarios' wavelength |
---|
563 | if (reff_driven_IRtoVIS_scenario .and. |
---|
564 | & (dustrefir.ne.9.3E-6)) then |
---|
565 | print*,'if reff_driven_IRtoVIS_scenario is used, then '// |
---|
566 | & 'the GCM IR reference wavelength should be the one '// |
---|
567 | & 'of the scenarios (dustiropacity=tes)' |
---|
568 | call abort_physic(modname, |
---|
569 | & "reff_driven_IRtoVIS_scenario requires tes wavelength",1) |
---|
570 | endif |
---|
571 | |
---|
572 | ! callddevil |
---|
573 | write(*,*)" dust lifted by dust devils ?" |
---|
574 | callddevil=.false. !default value |
---|
575 | call getin_p("callddevil",callddevil) |
---|
576 | write(*,*)" callddevil = ",callddevil |
---|
577 | |
---|
578 | ! Test of incompatibility: |
---|
579 | ! if dustdevil is used, then dustbin should be > 0 |
---|
580 | if (callddevil.and.(dustbin.lt.1)) then |
---|
581 | print*,'if dustdevil is used, then dustbin should > 0' |
---|
582 | call abort_physic(modname, |
---|
583 | & "callddevil requires dustbin > 0",1) |
---|
584 | endif |
---|
585 | |
---|
586 | ! sedimentation |
---|
587 | write(*,*) "Gravitationnal sedimentation ?" |
---|
588 | sedimentation=.true. ! default value |
---|
589 | call getin_p("sedimentation",sedimentation) |
---|
590 | write(*,*) " sedimentation = ",sedimentation |
---|
591 | ! activice |
---|
592 | write(*,*) "Radiatively active transported atmospheric ", |
---|
593 | & "water ice ?" |
---|
594 | activice=.false. ! default value |
---|
595 | call getin_p("activice",activice) |
---|
596 | write(*,*) " activice = ",activice |
---|
597 | ! water |
---|
598 | write(*,*) "Compute water cycle ?" |
---|
599 | water=.false. ! default value |
---|
600 | call getin_p("water",water) |
---|
601 | write(*,*) " water = ",water |
---|
602 | ! hdo |
---|
603 | write(*,*) "Compute hdo cycle ?" |
---|
604 | hdo=.false. ! default value |
---|
605 | call getin_p("hdo",hdo) |
---|
606 | write(*,*) " hdo = ",hdo |
---|
607 | |
---|
608 | write(*,*) "Use fractionation for hdo?" |
---|
609 | hdofrac=.true. ! default value |
---|
610 | call getin_p("hdofrac",hdofrac) |
---|
611 | write(*,*) " hdofrac = ",hdofrac |
---|
612 | |
---|
613 | ! Activeco2ice |
---|
614 | write(*,*) "Radiatively active transported atmospheric ", |
---|
615 | & "Co2 ice ?" |
---|
616 | activeco2ice=.false. ! default value |
---|
617 | call getin_p("activeco2ice",activeco2ice) |
---|
618 | write(*,*) " activeco2ice = ",activeco2ice |
---|
619 | |
---|
620 | ! sub-grid cloud fraction: fixed |
---|
621 | write(*,*) "Fixed cloud fraction?" |
---|
622 | CLFfixval=1.0 ! default value |
---|
623 | call getin_p("CLFfixval",CLFfixval) |
---|
624 | write(*,*) "CLFfixval=",CLFfixval |
---|
625 | ! sub-grid cloud fraction: varying |
---|
626 | write(*,*) "Use partial nebulosity?" |
---|
627 | CLFvarying=.false. ! default value |
---|
628 | call getin_p("CLFvarying",CLFvarying) |
---|
629 | write(*,*)"CLFvarying=",CLFvarying |
---|
630 | |
---|
631 | !CO2 clouds scheme? |
---|
632 | write(*,*) "Compute CO2 clouds (implies microphysical scheme)?" |
---|
633 | co2clouds=.false. ! default value |
---|
634 | call getin_p("co2clouds",co2clouds) |
---|
635 | write(*,*) " co2clouds = ",co2clouds |
---|
636 | !Can water ice particles serve as CCN for CO2clouds |
---|
637 | write(*,*) "Use water ice as CO2 clouds CCN ?" |
---|
638 | co2useh2o=.false. ! default value |
---|
639 | call getin_p("co2useh2o",co2useh2o) |
---|
640 | write(*,*) " co2useh2o = ",co2useh2o |
---|
641 | !Do we allow a supply of meteoritic paricles to serve as CO2 ice CCN? |
---|
642 | write(*,*) "Supply meteoritic particle for CO2 clouds ?" |
---|
643 | meteo_flux=.false. !Default value |
---|
644 | call getin_p("meteo_flux",meteo_flux) |
---|
645 | write(*,*) " meteo_flux = ",meteo_flux |
---|
646 | !Do we allow a sub-grid temperature distribution for the CO2 microphysics |
---|
647 | write(*,*) "sub-grid temperature distribution for CO2 clouds?" |
---|
648 | CLFvaryingCO2=.false. !Default value |
---|
649 | call getin_p("CLFvaryingCO2",CLFvaryingCO2) |
---|
650 | write(*,*) " CLFvaryingCO2 = ",CLFvaryingCO2 |
---|
651 | !Amplitude of the sub-grid temperature distribution for the CO2 microphysics |
---|
652 | write(*,*) "sub-grid temperature amplitude for CO2 clouds?" |
---|
653 | spantCO2=0 !Default value |
---|
654 | call getin_p("spantCO2",spantCO2) |
---|
655 | write(*,*) " spantCO2 = ",spantCO2 |
---|
656 | !Do you want to filter the sub-grid T distribution by a Saturation index? |
---|
657 | write(*,*) "filter sub-grid temperature by Saturation index?" |
---|
658 | satindexco2=.true. |
---|
659 | call getin_p("satindexco2",satindexco2) |
---|
660 | write(*,*) " satindexco2 = ",satindexco2 |
---|
661 | |
---|
662 | |
---|
663 | ! thermal inertia feedback |
---|
664 | write(*,*) "Activate the thermal inertia feedback ?" |
---|
665 | tifeedback=.false. ! default value |
---|
666 | call getin_p("tifeedback",tifeedback) |
---|
667 | write(*,*) " tifeedback = ",tifeedback |
---|
668 | |
---|
669 | ! Test of incompatibility: |
---|
670 | |
---|
671 | if (tifeedback.and..not.water) then |
---|
672 | print*,'if tifeedback is used,' |
---|
673 | print*,'water should be used too' |
---|
674 | call abort_physic(modname, |
---|
675 | & "tifeedback requires water",1) |
---|
676 | endif |
---|
677 | |
---|
678 | if (tifeedback.and..not.callsoil) then |
---|
679 | print*,'if tifeedback is used,' |
---|
680 | print*,'callsoil should be used too' |
---|
681 | call abort_physic(modname, |
---|
682 | & "tifeedback requires callsoil",1) |
---|
683 | endif |
---|
684 | |
---|
685 | if (activice.and..not.water) then |
---|
686 | print*,'if activice is used, water should be used too' |
---|
687 | call abort_physic(modname, |
---|
688 | & "activeice requires water",1) |
---|
689 | endif |
---|
690 | |
---|
691 | if (hdo.and..not.water) then |
---|
692 | print*,'if hdo is used, water should be used too' |
---|
693 | call abort_physic(modname, |
---|
694 | & "hd2 requires tracer",1) |
---|
695 | endif |
---|
696 | |
---|
697 | |
---|
698 | if (activeco2ice.and..not.co2clouds) then |
---|
699 | print*,'if activeco2ice is used, co2clouds should be used too' |
---|
700 | call abort_physic(modname, |
---|
701 | & "activeco2ice requires co2clouds",1) |
---|
702 | endif |
---|
703 | |
---|
704 | ! water ice clouds effective variance distribution for sedimentaion |
---|
705 | write(*,*) "Sed effective variance for water ice clouds ?" |
---|
706 | nuice_sed=0.45 |
---|
707 | call getin_p("nuice_sed",nuice_sed) |
---|
708 | write(*,*) "water_param nueff Sedimentation:", nuice_sed |
---|
709 | |
---|
710 | write(*,*) "Sed effective variance for CO2 clouds ?" |
---|
711 | nuiceco2_sed=0.45 |
---|
712 | call getin_p("nuiceco2_sed",nuiceco2_sed) |
---|
713 | write(*,*) "CO2 nueff Sedimentation:", nuiceco2_sed |
---|
714 | |
---|
715 | write(*,*) "REF effective variance for CO2 clouds ?" |
---|
716 | nuiceco2_ref=0.45 |
---|
717 | call getin_p("nuiceco2_ref",nuiceco2_ref) |
---|
718 | write(*,*) "CO2 nueff Sedimentation:", nuiceco2_ref |
---|
719 | |
---|
720 | write(*,*) "REF effective variance for water clouds ?" |
---|
721 | nuice_ref=0.1 |
---|
722 | call getin_p("nuice_ref",nuice_ref) |
---|
723 | write(*,*) "H2O nueff Sedimentation:", nuice_ref |
---|
724 | |
---|
725 | |
---|
726 | ! ccn factor if no scavenging |
---|
727 | write(*,*) "water param CCN reduc. factor ?" |
---|
728 | ccn_factor = 4.5 |
---|
729 | call getin_p("ccn_factor",ccn_factor) |
---|
730 | write(*,*)" ccn_factor = ",ccn_factor |
---|
731 | write(*,*)"Careful: only used when microphys=F, otherwise" |
---|
732 | write(*,*)"the contact parameter is used instead;" |
---|
733 | |
---|
734 | ! microphys |
---|
735 | write(*,*)"Microphysical scheme for water-ice clouds?" |
---|
736 | microphys=.false. ! default value |
---|
737 | call getin_p("microphys",microphys) |
---|
738 | write(*,*)" microphys = ",microphys |
---|
739 | |
---|
740 | ! supersat |
---|
741 | write(*,*)"Allow super-saturation of water vapor?" |
---|
742 | supersat=.true. ! default value |
---|
743 | call getin_p("supersat",supersat) |
---|
744 | write(*,*)"supersat = ",supersat |
---|
745 | |
---|
746 | ! microphysical parameter contact |
---|
747 | write(*,*) "water contact parameter ?" |
---|
748 | mteta = 0.95 ! default value |
---|
749 | temp_dependent_m = .false. ! default value |
---|
750 | call getin_p("temp_dependent_m",temp_dependent_m) |
---|
751 | if (temp_dependent_m) then !(JN 2023) |
---|
752 | print*,'You have chosen a temperature-dependent water' |
---|
753 | print*,'contact parameter ! From Maattanen et al. 2014' |
---|
754 | else if (.not.temp_dependent_m) then |
---|
755 | print*,'Water contact parameter is constant' |
---|
756 | call getin_p("mteta",mteta) |
---|
757 | write(*,*) "mteta = ", mteta |
---|
758 | endif |
---|
759 | |
---|
760 | ! Adaptative timestep for cloud microphysics (JN 2023) |
---|
761 | write(*,*)"Adaptative timestep for cloud", |
---|
762 | & " microphysics ? (default is false)" |
---|
763 | cloud_adapt_ts=.false. ! default value |
---|
764 | call getin_p("cloud_adapt_ts",cloud_adapt_ts) |
---|
765 | write(*,*)"cloud_adapt_ts= ",cloud_adapt_ts |
---|
766 | |
---|
767 | ! Test of incompatibility: |
---|
768 | ! If you want the adaptative timestep, you should use the |
---|
769 | ! temperature dependent contact parameter (otherwise the |
---|
770 | ! global water cycle is nuts !) |
---|
771 | ! However one can use the temperature dependent contact parameter |
---|
772 | ! without the adaptative timestep (MCD6.1 configuration) |
---|
773 | |
---|
774 | if (cloud_adapt_ts.and..not.temp_dependent_m) then |
---|
775 | print*,'Water cycle v6 : if cloud_adapt_ts is used' |
---|
776 | print*,'then temp_dependent_m must be used!' |
---|
777 | print*,'Otherwise the water cycle is unrealistic.' |
---|
778 | call abort_physic(modname, |
---|
779 | & "cloud_adapt_ts requires temp_dependent_m",1) |
---|
780 | endif |
---|
781 | ! scavenging |
---|
782 | write(*,*)"Dust scavenging by H2O/CO2 snowfall ?" |
---|
783 | scavenging=.false. ! default value |
---|
784 | call getin_p("scavenging",scavenging) |
---|
785 | write(*,*)" scavenging = ",scavenging |
---|
786 | |
---|
787 | |
---|
788 | ! Test of incompatibility: |
---|
789 | ! if scavenging is used, then dustbin should be > 0 |
---|
790 | |
---|
791 | if ((microphys.and..not.doubleq).or. |
---|
792 | & (microphys.and..not.water)) then |
---|
793 | print*,'if microphys is used, then doubleq,' |
---|
794 | print*,'and water must be used!' |
---|
795 | call abort_physic(modname, |
---|
796 | & "microphys requires water and doubleq",1) |
---|
797 | endif |
---|
798 | if (microphys.and..not.scavenging) then |
---|
799 | print*,'' |
---|
800 | print*,'----------------WARNING-----------------' |
---|
801 | print*,'microphys is used without scavenging !!!' |
---|
802 | print*,'----------------WARNING-----------------' |
---|
803 | print*,'' |
---|
804 | endif |
---|
805 | |
---|
806 | if ((scavenging.and..not.microphys).or. |
---|
807 | & (scavenging.and.(dustbin.lt.1)))then |
---|
808 | print*,'if scavenging is used, then microphys' |
---|
809 | print*,'must be used!' |
---|
810 | call abort_physic(modname, |
---|
811 | & "scavenging requires microphys",1) |
---|
812 | endif |
---|
813 | |
---|
814 | ! Instantaneous scavenging by CO2 |
---|
815 | ! -> expected to be replaced by scavenging with microphysics (flag scavenging) one day |
---|
816 | write(*,*)"Dust scavenging by instantaneous CO2 snowfall ?" |
---|
817 | scavco2cond=.false. ! default value |
---|
818 | call getin_p("scavco2cond",scavco2cond) |
---|
819 | write(*,*)" scavco2cond = ",scavco2cond |
---|
820 | ! Test of incompatibility: |
---|
821 | ! if scavco2cond is used, then dustbin should be > 0 |
---|
822 | if (scavco2cond.and.(dustbin.lt.1))then |
---|
823 | print*,'if scavco2cond is used, then dustbin should be > 0' |
---|
824 | call abort_physic(modname, |
---|
825 | & "scavco2cond requires dustbin > 0",1) |
---|
826 | endif |
---|
827 | ! if co2clouds is used, then there is no need for scavco2cond |
---|
828 | if (co2clouds.and.scavco2cond) then |
---|
829 | print*,'' |
---|
830 | print*,'----------------WARNING-----------------' |
---|
831 | print*,' microphys scavenging is used so ' |
---|
832 | print*,' no need for scavco2cond !!! ' |
---|
833 | print*,'----------------WARNING-----------------' |
---|
834 | print*,'' |
---|
835 | call abort_physic(modname, |
---|
836 | & "incompatible co2cloud and scavco2cond options",1) |
---|
837 | endif |
---|
838 | |
---|
839 | ! Test of incompatibility: |
---|
840 | |
---|
841 | write(*,*) "Permanent water caps at poles ?", |
---|
842 | & " .true. is RECOMMENDED" |
---|
843 | write(*,*) "(with .true., North cap is a source of water ", |
---|
844 | & "and South pole is a cold trap)" |
---|
845 | caps=.true. ! default value |
---|
846 | call getin_p("caps",caps) |
---|
847 | write(*,*) " caps = ",caps |
---|
848 | |
---|
849 | ! JN : now separated between albedo_h2o_cap and |
---|
850 | ! albedo_h2o_frost. Retrocompatible with old |
---|
851 | ! callphys.def with albedo_h2o_ice |
---|
852 | write(*,*) "water ice albedo ? Old settings use ", |
---|
853 | & "albedo_h2o_ice, new settings use ", |
---|
854 | & "albedo_h2o_cap and albedo_h2o_frost " |
---|
855 | albedo_h2o_cap=0.35 |
---|
856 | albedo_h2o_frost=0.35 |
---|
857 | call getin_p("albedo_h2o_ice",albedo_h2o_cap) |
---|
858 | albedo_h2o_frost=albedo_h2o_cap |
---|
859 | call getin_p("albedo_h2o_cap",albedo_h2o_cap) |
---|
860 | write(*,*) " albedo_h2o_cap = ",albedo_h2o_cap |
---|
861 | call getin_p("albedo_h2o_frost",albedo_h2o_frost) |
---|
862 | write(*,*) " albedo_h2o_frost = ",albedo_h2o_frost |
---|
863 | |
---|
864 | ! Northern polar cap albedo (JN 2021) |
---|
865 | write(*,*)"Watercaptag albedo is unchanged by water frost", |
---|
866 | & " deposition (default is false)" |
---|
867 | cst_cap_albedo=.false. ! default value |
---|
868 | call getin_p("cst_cap_albedo",cst_cap_albedo) |
---|
869 | write(*,*)"cst_cap_albedo = ",cst_cap_albedo |
---|
870 | |
---|
871 | ! Watercap evolution & refill (with discriminated albedos) (JN 2021) |
---|
872 | write(*,*)"Watercap is replenished by water frost", |
---|
873 | & " accumulation (default is false)" |
---|
874 | refill_watercap=.false. ! default value |
---|
875 | call getin_p("refill_watercap",refill_watercap) |
---|
876 | write(*,*)"refill_watercap= ",refill_watercap |
---|
877 | ! frost thickness threshold for refill_watercap (ice metamorphism) |
---|
878 | write(*,*) "frost thickness threshold for metamorphism ?", |
---|
879 | & "ie converted into watercap", |
---|
880 | & "only if refill_watercap is .true." |
---|
881 | frost_metam_threshold=0.05 ! (i.e 0.05 kg.m-2) |
---|
882 | call getin_p("frost_metam_threshold", |
---|
883 | & frost_metam_threshold) |
---|
884 | write(*,*) " frost_metam_threshold = ", |
---|
885 | & frost_metam_threshold |
---|
886 | |
---|
887 | |
---|
888 | ! inert_h2o_ice |
---|
889 | write(*,*) "water ice thermal inertia ?" |
---|
890 | inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2) |
---|
891 | call getin_p("inert_h2o_ice",inert_h2o_ice) |
---|
892 | write(*,*) " inert_h2o_ice = ",inert_h2o_ice |
---|
893 | ! frost_albedo_threshold |
---|
894 | write(*,*) "frost thickness threshold for albedo ?" |
---|
895 | frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2) |
---|
896 | call getin_p("frost_albedo_threshold", |
---|
897 | & frost_albedo_threshold) |
---|
898 | write(*,*) " frost_albedo_threshold = ", |
---|
899 | & frost_albedo_threshold |
---|
900 | |
---|
901 | ! TMP: old_wsublimation_scheme |
---|
902 | write(*,*) "Old water sublimation scheme?" |
---|
903 | old_wsublimation_scheme = .true. |
---|
904 | call getin_p("old_wsublimation_scheme",old_wsublimation_scheme) |
---|
905 | write(*,*) "old_wsublimation_scheme",old_wsublimation_scheme |
---|
906 | |
---|
907 | ! call Titus crocus line -- DEFAULT IS NONE |
---|
908 | write(*,*) "Titus crocus line ?" |
---|
909 | tituscap=.false. ! default value |
---|
910 | call getin_p("tituscap",tituscap) |
---|
911 | write(*,*) "tituscap",tituscap |
---|
912 | |
---|
913 | ! Chemistry: |
---|
914 | write(*,*) "photochemistry: include chemical species" |
---|
915 | photochem=.false. ! default value |
---|
916 | call getin_p("photochem",photochem) |
---|
917 | write(*,*) " photochem = ",photochem |
---|
918 | |
---|
919 | write(*,*) "Compute chemistry (if photochem is .true.)", |
---|
920 | & "every ichemistry physics step (default: ichemistry=1)" |
---|
921 | ichemistry=1 |
---|
922 | call getin_p("ichemistry",ichemistry) |
---|
923 | write(*,*) " ichemistry = ",ichemistry |
---|
924 | |
---|
925 | |
---|
926 | ! SCATTERERS |
---|
927 | write(*,*) "how many scatterers?" |
---|
928 | naerkind=1 ! default value |
---|
929 | call getin_p("naerkind",naerkind) |
---|
930 | write(*,*)" naerkind = ",naerkind |
---|
931 | |
---|
932 | ! Test of incompatibility |
---|
933 | c Logical tests for radiatively active water-ice clouds: |
---|
934 | IF ( (activice.AND.(.NOT.water)).OR. |
---|
935 | & (activice.AND.(naerkind.LT.2)) ) THEN |
---|
936 | WRITE(*,*) 'If activice is TRUE, water has to be set' |
---|
937 | WRITE(*,*) 'to TRUE, and "naerkind" must be at least' |
---|
938 | WRITE(*,*) 'equal to 2.' |
---|
939 | call abort_physic(modname, |
---|
940 | & "radiatively active dust and water"// |
---|
941 | & " require naerkind > 1",1) |
---|
942 | ENDIF |
---|
943 | |
---|
944 | !------------------------------------------ |
---|
945 | !------------------------------------------ |
---|
946 | ! once naerkind is known allocate arrays |
---|
947 | ! -- we do it here and not in phys_var_init |
---|
948 | ! -- because we need to know naerkind first |
---|
949 | CALL ini_scatterers(ngrid,nlayer) |
---|
950 | !------------------------------------------ |
---|
951 | !------------------------------------------ |
---|
952 | |
---|
953 | |
---|
954 | c Please name the different scatterers here ---------------- |
---|
955 | name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile |
---|
956 | IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme |
---|
957 | |
---|
958 | if (nq.gt.1) then |
---|
959 | ! trick to avoid problems compiling with 1 tracer |
---|
960 | ! and picky compilers who know name_iaer(2) is out of bounds |
---|
961 | j=2 |
---|
962 | IF (rdstorm.AND..NOT.activice.AND..NOT.topflows) then |
---|
963 | name_iaer(j) = "stormdust_doubleq" !! storm dust two-moment scheme |
---|
964 | j = j+1 |
---|
965 | END IF |
---|
966 | |
---|
967 | IF (rdstorm.AND.water.AND.activice.AND..NOT.topflows) then |
---|
968 | name_iaer(j) = "stormdust_doubleq" |
---|
969 | j = j+1 |
---|
970 | END IF |
---|
971 | |
---|
972 | IF (topflows.AND..NOT.activice.AND..NOT.rdstorm) then |
---|
973 | name_iaer(j) = "topdust_doubleq" !! storm dust two-moment scheme |
---|
974 | j = j+1 |
---|
975 | END IF |
---|
976 | |
---|
977 | IF (topflows.AND.water.AND.activice.AND..NOT.rdstorm) then |
---|
978 | name_iaer(j) = "topdust_doubleq" |
---|
979 | j = j+1 |
---|
980 | END IF |
---|
981 | |
---|
982 | IF (rdstorm.AND.topflows.AND..NOT.activice) THEN |
---|
983 | name_iaer(j) = "stormdust_doubleq" |
---|
984 | name_iaer(j+1) = "topdust_doubleq" |
---|
985 | j = j+2 |
---|
986 | ENDIF |
---|
987 | |
---|
988 | IF (rdstorm.AND.topflows.AND.water.AND.activice) THEN |
---|
989 | name_iaer(j) = "stormdust_doubleq" |
---|
990 | name_iaer(j+1) = "topdust_doubleq" |
---|
991 | j = j+2 |
---|
992 | ENDIF |
---|
993 | |
---|
994 | IF (water.AND.activice) then |
---|
995 | name_iaer(j) = "h2o_ice" !! radiatively-active clouds |
---|
996 | j = j+1 |
---|
997 | END IF |
---|
998 | |
---|
999 | IF (co2clouds.AND.activeco2ice) then |
---|
1000 | name_iaer(j) = "co2_ice" !! radiatively-active co2 clouds |
---|
1001 | j = j+1 |
---|
1002 | ENDIF |
---|
1003 | |
---|
1004 | IF (submicron.AND.active) then |
---|
1005 | name_iaer(j) = "dust_submicron" !! JBM experimental stuff |
---|
1006 | j = j+1 |
---|
1007 | ENDIF |
---|
1008 | endif ! of if (nq.gt.1) |
---|
1009 | c ---------------------------------------------------------- |
---|
1010 | |
---|
1011 | ! Adsorption |
---|
1012 | adsorption_soil = .false. ! default value |
---|
1013 | call getin_p("adsorption_soil",adsorption_soil) |
---|
1014 | if (adsorption_soil .and. (.not. water)) then |
---|
1015 | write(*,*)"Adsorption can be run only if water = True" |
---|
1016 | call abort_physic(modname, |
---|
1017 | & "Adsorption must be used with water = true",1) |
---|
1018 | endif |
---|
1019 | |
---|
1020 | choice_ads = 0 ! default value (no adsorption) |
---|
1021 | if (adsorption_soil) choice_ads = 1 |
---|
1022 | call getin_p("choice_ads",choice_ads) |
---|
1023 | |
---|
1024 | c ---------------------------------------------------------- |
---|
1025 | |
---|
1026 | ! THERMOSPHERE |
---|
1027 | |
---|
1028 | write(*,*) "call thermosphere ?" |
---|
1029 | callthermos=.false. ! default value |
---|
1030 | call getin_p("callthermos",callthermos) |
---|
1031 | write(*,*) " callthermos = ",callthermos |
---|
1032 | |
---|
1033 | |
---|
1034 | write(*,*) " water included without cycle ", |
---|
1035 | & "(only if water=.false.)" |
---|
1036 | thermoswater=.false. ! default value |
---|
1037 | call getin_p("thermoswater",thermoswater) |
---|
1038 | write(*,*) " thermoswater = ",thermoswater |
---|
1039 | |
---|
1040 | write(*,*) "call thermal conduction ?", |
---|
1041 | & " (only if callthermos=.true.)" |
---|
1042 | callconduct=.false. ! default value |
---|
1043 | call getin_p("callconduct",callconduct) |
---|
1044 | write(*,*) " callconduct = ",callconduct |
---|
1045 | |
---|
1046 | write(*,*) "call EUV heating ?", |
---|
1047 | & " (only if callthermos=.true.)" |
---|
1048 | calleuv=.false. ! default value |
---|
1049 | call getin_p("calleuv",calleuv) |
---|
1050 | write(*,*) " calleuv = ",calleuv |
---|
1051 | |
---|
1052 | write(*,*) "call molecular viscosity ?", |
---|
1053 | & " (only if callthermos=.true.)" |
---|
1054 | callmolvis=.false. ! default value |
---|
1055 | call getin_p("callmolvis",callmolvis) |
---|
1056 | write(*,*) " callmolvis = ",callmolvis |
---|
1057 | |
---|
1058 | write(*,*) "call molecular diffusion ?", |
---|
1059 | & " (only if callthermos=.true.)" |
---|
1060 | callmoldiff=.false. ! default value |
---|
1061 | call getin_p("callmoldiff",callmoldiff) |
---|
1062 | write(*,*) " callmoldiff = ",callmoldiff |
---|
1063 | |
---|
1064 | |
---|
1065 | write(*,*) "call thermospheric photochemistry ?", |
---|
1066 | & " (only if callthermos=.true.)" |
---|
1067 | thermochem=.false. ! default value |
---|
1068 | call getin_p("thermochem",thermochem) |
---|
1069 | write(*,*) " thermochem = ",thermochem |
---|
1070 | |
---|
1071 | write(*,*) "Method to include solar variability" |
---|
1072 | write(*,*) "0-> fixed value of E10.7 (fixed_euv_value); ", |
---|
1073 | & "1-> daily evolution of E10.7 (for given solvaryear)" |
---|
1074 | solvarmod=1 |
---|
1075 | call getin_p("solvarmod",solvarmod) |
---|
1076 | write(*,*) " solvarmod = ",solvarmod |
---|
1077 | |
---|
1078 | write(*,*) "Fixed euv (for solvarmod==0) 10.7 value?" |
---|
1079 | write(*,*) " (min=80 , ave=140, max=320)" |
---|
1080 | fixed_euv_value=140 ! default value |
---|
1081 | call getin_p("fixed_euv_value",fixed_euv_value) |
---|
1082 | write(*,*) " fixed_euv_value = ",fixed_euv_value |
---|
1083 | |
---|
1084 | write(*,*) "Solar variability as observed for MY: " |
---|
1085 | write(*,*) "Only if solvarmod=1" |
---|
1086 | solvaryear=24 |
---|
1087 | call getin_p("solvaryear",solvaryear) |
---|
1088 | write(*,*) " solvaryear = ",solvaryear |
---|
1089 | |
---|
1090 | write(*,*) "UV heating efficiency:", |
---|
1091 | & "measured values between 0.19 and 0.23 (Fox et al. 1996)", |
---|
1092 | & "lower values may be used to compensate low 15 um cooling" |
---|
1093 | euveff=0.21 !default value |
---|
1094 | call getin_p("euveff",euveff) |
---|
1095 | write(*,*) " euveff = ", euveff |
---|
1096 | |
---|
1097 | |
---|
1098 | if (.not.callthermos) then |
---|
1099 | if (thermoswater) then |
---|
1100 | print*,'if thermoswater is set, callthermos must be true' |
---|
1101 | call abort_physic(modname, |
---|
1102 | & "thermoswater requires callthermos",1) |
---|
1103 | endif |
---|
1104 | if (callconduct) then |
---|
1105 | print*,'if callconduct is set, callthermos must be true' |
---|
1106 | call abort_physic(modname, |
---|
1107 | & "callconduct requires callthermos",1) |
---|
1108 | endif |
---|
1109 | if (calleuv) then |
---|
1110 | print*,'if calleuv is set, callthermos must be true' |
---|
1111 | call abort_physic(modname, |
---|
1112 | & "calleuv requires callthermos",1) |
---|
1113 | endif |
---|
1114 | if (callmolvis) then |
---|
1115 | print*,'if callmolvis is set, callthermos must be true' |
---|
1116 | call abort_physic(modname, |
---|
1117 | & "callmolvis requires callthermos",1) |
---|
1118 | endif |
---|
1119 | if (callmoldiff) then |
---|
1120 | print*,'if callmoldiff is set, callthermos must be true' |
---|
1121 | call abort_physic(modname, |
---|
1122 | & "callmoldiff requires callthermos",1) |
---|
1123 | endif |
---|
1124 | if (thermochem) then |
---|
1125 | print*,'if thermochem is set, callthermos must be true' |
---|
1126 | call abort_physic(modname, |
---|
1127 | & "thermochem requires callthermos",1) |
---|
1128 | endif |
---|
1129 | endif |
---|
1130 | |
---|
1131 | ! Test of incompatibility: |
---|
1132 | ! if photochem is used, then water should be used too |
---|
1133 | |
---|
1134 | if (photochem.and..not.water) then |
---|
1135 | print*,'if photochem is used, water should be used too' |
---|
1136 | call abort_physic(modname, |
---|
1137 | & "photochem requires water",1) |
---|
1138 | endif |
---|
1139 | |
---|
1140 | ! if callthermos is used, then thermoswater should be used too |
---|
1141 | ! (if water not used already) |
---|
1142 | |
---|
1143 | if (callthermos .and. .not.water) then |
---|
1144 | if (callthermos .and. .not.thermoswater) then |
---|
1145 | print*,'if callthermos is used, water or thermoswater |
---|
1146 | & should be used too' |
---|
1147 | call abort_physic(modname, |
---|
1148 | & "callthermos requires water or thermoswater",1) |
---|
1149 | endif |
---|
1150 | endif |
---|
1151 | |
---|
1152 | PRINT*,'--------------------------------------------' |
---|
1153 | PRINT* |
---|
1154 | PRINT* |
---|
1155 | ELSE |
---|
1156 | write(*,*) |
---|
1157 | write(*,*) 'Cannot read file callphys.def. Is it here ?' |
---|
1158 | call abort_physic(modname, |
---|
1159 | & "missing callphys.def file",1) |
---|
1160 | ENDIF |
---|
1161 | |
---|
1162 | 8000 FORMAT(t5,a12,l8) |
---|
1163 | 8001 FORMAT(t5,a12,i8) |
---|
1164 | |
---|
1165 | PRINT* |
---|
1166 | PRINT*,'conf_phys: daysec',daysec |
---|
1167 | PRINT* |
---|
1168 | PRINT*,'conf_phys: The radiative transfer is computed:' |
---|
1169 | PRINT*,' each ',iradia,' physical time-step' |
---|
1170 | PRINT*,' or each ',iradia*dtphys,' seconds' |
---|
1171 | PRINT* |
---|
1172 | ! -------------------------------------------------------------- |
---|
1173 | ! Managing the Longwave radiative transfer |
---|
1174 | ! -------------------------------------------------------------- |
---|
1175 | |
---|
1176 | ! In most cases, the run just use the following values : |
---|
1177 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
1178 | callemis=.true. |
---|
1179 | ! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01 |
---|
1180 | ilwd=1 |
---|
1181 | ilwn=1 !2 |
---|
1182 | ilwb=1 !2 |
---|
1183 | linear=.true. |
---|
1184 | ncouche=3 |
---|
1185 | alphan=0.4 |
---|
1186 | semi=0 |
---|
1187 | |
---|
1188 | ! BUT people working hard on the LW may want to read them in 'radia.def' |
---|
1189 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
1190 | !$OMP MASTER |
---|
1191 | OPEN(99,file='radia.def',status='old',form='formatted' |
---|
1192 | . ,iostat=ierr) |
---|
1193 | IF(ierr.EQ.0) THEN |
---|
1194 | write(*,*) 'conf_phys: Reading radia.def !!!' |
---|
1195 | READ(99,fmt='(a)') ch1 |
---|
1196 | READ(99,*) callemis |
---|
1197 | WRITE(*,8000) ch1,callemis |
---|
1198 | |
---|
1199 | READ(99,fmt='(a)') ch1 |
---|
1200 | READ(99,*) iradia |
---|
1201 | WRITE(*,8001) ch1,iradia |
---|
1202 | |
---|
1203 | READ(99,fmt='(a)') ch1 |
---|
1204 | READ(99,*) ilwd |
---|
1205 | WRITE(*,8001) ch1,ilwd |
---|
1206 | |
---|
1207 | READ(99,fmt='(a)') ch1 |
---|
1208 | READ(99,*) ilwn |
---|
1209 | WRITE(*,8001) ch1,ilwn |
---|
1210 | |
---|
1211 | READ(99,fmt='(a)') ch1 |
---|
1212 | READ(99,*) linear |
---|
1213 | WRITE(*,8000) ch1,linear |
---|
1214 | |
---|
1215 | READ(99,fmt='(a)') ch1 |
---|
1216 | READ(99,*) ncouche |
---|
1217 | WRITE(*,8001) ch1,ncouche |
---|
1218 | |
---|
1219 | READ(99,fmt='(a)') ch1 |
---|
1220 | READ(99,*) alphan |
---|
1221 | WRITE(*,*) ch1,alphan |
---|
1222 | |
---|
1223 | READ(99,fmt='(a)') ch1 |
---|
1224 | READ(99,*) ilwb |
---|
1225 | WRITE(*,8001) ch1,ilwb |
---|
1226 | |
---|
1227 | |
---|
1228 | READ(99,fmt='(a)') ch1 |
---|
1229 | READ(99,'(l1)') callg2d |
---|
1230 | WRITE(*,8000) ch1,callg2d |
---|
1231 | |
---|
1232 | READ(99,fmt='(a)') ch1 |
---|
1233 | READ(99,*) semi |
---|
1234 | WRITE(*,*) ch1,semi |
---|
1235 | end if |
---|
1236 | CLOSE(99) |
---|
1237 | !$OMP END MASTER |
---|
1238 | call bcast(ch1) |
---|
1239 | call bcast(callemis) |
---|
1240 | call bcast(iradia) |
---|
1241 | call bcast(ilwd) |
---|
1242 | call bcast(ilwn) |
---|
1243 | call bcast(linear) |
---|
1244 | call bcast(ncouche) |
---|
1245 | call bcast(alphan) |
---|
1246 | call bcast(ilwb) |
---|
1247 | call bcast(callg2d) |
---|
1248 | call bcast(semi) |
---|
1249 | |
---|
1250 | END SUBROUTINE conf_phys |
---|
1251 | |
---|
1252 | END MODULE conf_phys_mod |
---|