source: trunk/LMDZ.MARS/libf/phymars/conf_phys.F @ 3712

Last change on this file since 3712 was 3712, checked in by emillour, 3 months ago

Mars PCM:
Add extra checking in nirco2abs to ensure having nircorr==1 when o and
co2 tracers are available.
EM

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