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

Last change on this file since 1381 was 1381, checked in by emillour, 10 years ago

Mars GCM:

EM

File size: 27.7 KB
Line 
1      SUBROUTINE conf_phys(ngrid,nlayer,nq)
2 
3!=======================================================================
4!
5!   purpose:
6!   -------
7!
8!   Initialisation for the physical parametrisations of the LMD
9!   martian atmospheric general circulation modele.
10!
11!   author: Frederic Hourdin 15 / 10 /93
12!   -------
13!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
14!             Ehouarn Millour (oct. 2008) tracers are now identified
15!              by their names and may not be contiguously
16!              stored in the q(:,:,:,:) array
17!             E.M. (june 2009) use getin routine to load parameters
18!             adapted to the mesoscale use - Aymeric Spiga - 01/2007-07/2011
19!             separated inifis into conf_phys and phys_state_var_init (A. Spiga)
20!
21!
22!   arguments:
23!   ----------
24!
25!   input:
26!   ------
27!
28!    nq                    Number of tracers
29!
30!=======================================================================
31!
32!-----------------------------------------------------------------------
33!   declarations:
34!   -------------
35! to use  'getin'
36      USE ioipsl_getincom, only : getin
37      use tracer_mod, only : nuice_sed, ccn_factor
38      use surfdat_h, only: albedo_h2o_ice, inert_h2o_ice,
39     &                     frost_albedo_threshold
40      use control_mod, only: ecritphy
41      use planete_h
42      USE comcstfi_h, only: daysec,dtphys
43      use dimradmars_mod, only: naerkind, name_iaer,
44     &                      ini_scatterers,tauvis
45
46      IMPLICIT NONE
47#include "callkeys.h"
48#include "datafile.h"
49#include "microphys.h"
50
51      INTEGER,INTENT(IN) :: ngrid,nlayer,nq
52      INTEGER ig,ierr,j
53 
54      CHARACTER ch1*12
55      CHARACTER ch80*80
56
57      ! read in 'ecritphy' (frequency of calls to physics, in dynamical steps)
58      ! (also done in dyn3d/defrun_new but not in LMDZ.COMMON)
59      call getin("ecritphy",ecritphy)
60     
61! --------------------------------------------------------------
62!  Reading the "callphys.def" file controlling some key options
63! --------------------------------------------------------------
64     
65      ! check that 'callphys.def' file is around
66      OPEN(99,file='callphys.def',status='old',form='formatted'
67     &     ,iostat=ierr)
68      CLOSE(99)
69     
70      IF(ierr.EQ.0) THEN
71         PRINT*
72         PRINT*
73         PRINT*,'--------------------------------------------'
74         PRINT*,' conf_phys: Parameters for the physics (callphys.def)'
75         PRINT*,'--------------------------------------------'
76
77         write(*,*) "Directory where external input files are:"
78         datafile="/u/lmdz/WWW/planets/mars/datadir"
79         call getin("datadir",datafile) ! default path
80         write(*,*) " datafile = ",trim(datafile)
81
82         write(*,*) "Run with or without tracer transport ?"
83         tracer=.false. ! default value
84         call getin("tracer",tracer)
85         write(*,*) " tracer = ",tracer
86
87         write(*,*) "Diurnal cycle ?"
88         write(*,*) "(if diurnal=False, diurnal averaged solar heating)"
89         diurnal=.true. ! default value
90         call getin("diurnal",diurnal)
91         write(*,*) " diurnal = ",diurnal
92
93         write(*,*) "Seasonal cycle ?"
94         write(*,*) "(if season=False, Ls stays constant, to value ",
95     &   "set in 'start'"
96         season=.true. ! default value
97         call getin("season",season)
98         write(*,*) " season = ",season
99
100         write(*,*) "Write some extra output to the screen ?"
101         lwrite=.false. ! default value
102         call getin("lwrite",lwrite)
103         write(*,*) " lwrite = ",lwrite
104
105         write(*,*) "Save statistics in file stats.nc ?"
106#ifdef MESOSCALE
107         callstats=.false. ! default value
108#else
109         callstats=.true. ! default value
110#endif
111         call getin("callstats",callstats)
112         write(*,*) " callstats = ",callstats
113
114         write(*,*) "Save EOF profiles in file 'profiles' for ",
115     &              "Climate Database?"
116         calleofdump=.false. ! default value
117         call getin("calleofdump",calleofdump)
118         write(*,*) " calleofdump = ",calleofdump
119
120         write(*,*) "Dust scenario: 1=constant dust (read from startfi",
121     &   " or set as tauvis); 2=Viking scenario; =3 MGS scenario,",
122     &   "=6 cold (low dust) scenario; =7 warm (high dust) scenario ",
123     &   "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation"
124         iaervar=3 ! default value
125         call getin("iaervar",iaervar)
126         write(*,*) " iaervar = ",iaervar
127
128         write(*,*) "Reference (visible) dust opacity at 610 Pa ",
129     &   "(matters only if iaervar=1)"
130         ! NB: default value of tauvis is set/read in startfi.nc file
131         call getin("tauvis",tauvis)
132         write(*,*) " tauvis = ",tauvis
133
134         write(*,*) "Dust vertical distribution:"
135         write(*,*) "(=1 top set by topdustref parameter;",
136     & " =2 Viking scenario; =3 MGS scenario)"
137         iddist=3 ! default value
138         call getin("iddist",iddist)
139         write(*,*) " iddist = ",iddist
140
141         write(*,*) "Dust top altitude (km). (Matters only if iddist=1)"
142         topdustref= 90.0 ! default value
143         call getin("topdustref",topdustref)
144         write(*,*) " topdustref = ",topdustref
145
146         write(*,*) "Prescribed surface thermal flux (H/(rho*cp),K m/s)"
147         tke_heat_flux=0. ! default value
148         call getin("tke_heat_flux",tke_heat_flux)
149         write(*,*) " tke_heat_flux = ",tke_heat_flux
150         write(*,*) " 0 means the usual schemes are computing"
151
152         write(*,*) "call radiative transfer ?"
153         callrad=.true. ! default value
154         call getin("callrad",callrad)
155         write(*,*) " callrad = ",callrad
156
157         write(*,*) "call slope insolation scheme ?",
158     &              "(matters only if callrad=T)"
159#ifdef MESOSCALE
160         callslope=.true. ! default value
161#else
162         callslope=.false. ! default value (not supported yet)
163#endif
164         call getin("callslope",callslope)
165         write(*,*) " callslope = ",callslope
166
167         write(*,*) "call NLTE radiative schemes ?",
168     &              "(matters only if callrad=T)"
169         callnlte=.false. ! default value
170         call getin("callnlte",callnlte)
171         write(*,*) " callnlte = ",callnlte
172         
173         nltemodel=0    !default value
174         write(*,*) "NLTE model?"
175         write(*,*) "0 -> old model, static O"
176         write(*,*) "1 -> old model, dynamic O"
177         write(*,*) "2 -> new model"
178         write(*,*) "(matters only if callnlte=T)"
179         call getin("nltemodel",nltemodel)
180         write(*,*) " nltemodel = ",nltemodel
181
182         write(*,*) "call CO2 NIR absorption ?",
183     &              "(matters only if callrad=T)"
184         callnirco2=.false. ! default value
185         call getin("callnirco2",callnirco2)
186         write(*,*) " callnirco2 = ",callnirco2
187
188         write(*,*) "New NIR NLTE correction ?",
189     $              "0-> old model (no correction)",
190     $              "1-> new correction",
191     $              "(matters only if callnirco2=T)"
192#ifdef MESOSCALE
193         nircorr=0      !default value. this is OK below 60 km.
194#else
195         nircorr=0      !default value
196#endif
197         call getin("nircorr",nircorr)
198         write(*,*) " nircorr = ",nircorr
199
200         write(*,*) "call turbulent vertical diffusion ?"
201         calldifv=.true. ! default value
202         call getin("calldifv",calldifv)
203         write(*,*) " calldifv = ",calldifv
204
205         write(*,*) "call thermals ?"
206         calltherm=.false. ! default value
207         call getin("calltherm",calltherm)
208         write(*,*) " calltherm = ",calltherm
209
210         write(*,*) "call convective adjustment ?"
211         calladj=.true. ! default value
212         call getin("calladj",calladj)
213         write(*,*) " calladj = ",calladj
214
215         if (calltherm .and. calladj) then
216          print*,'!!! PLEASE NOTE !!!'
217          print*,'convective adjustment is on'
218          print*,'but since thermal plume model is on'
219          print*,'convadj is only activated above the PBL'
220         endif
221       
222         write(*,*) "used latest version of yamada scheme?"
223         callyamada4=.true. ! default value
224         call getin("callyamada4",callyamada4)
225         write(*,*) " callyamada4 = ",callyamada4
226
227         if (calltherm .and. .not.callyamada4) then
228          print*,'!!!! WARNING WARNING WARNING !!!!'
229          print*,'if calltherm=T we strongly advise that '
230          print*,'you set the flag callyamada4 to T '
231          print*,'!!!! WARNING WARNING WARNING !!!!'
232         endif
233 
234         write(*,*) "call Richardson-based surface layer ?"
235         callrichsl=.false. ! default value
236         call getin("callrichsl",callrichsl)
237         write(*,*) " callrichsl = ",callrichsl
238
239         if (calltherm .and. .not.callrichsl) then
240          print*,'WARNING WARNING WARNING'
241          print*,'if calltherm=T we strongly advise that '
242          print*,'you use the new surface layer scheme '
243          print*,'by setting callrichsl=T '
244         endif
245
246         if (calladj .and. callrichsl .and. (.not. calltherm)) then
247          print*,'You should not be calling the convective adjustment
248     & scheme with the Richardson surface-layer and without the thermals
249     &. This approach is not
250     & physically consistent and can lead to unrealistic friction
251     & values.'
252          print*,'If you want to use the Ri. surface-layer, either
253     & activate thermals OR de-activate the convective adjustment.'
254          stop
255         endif
256
257         write(*,*) "call CO2 condensation ?"
258         callcond=.true. ! default value
259         call getin("callcond",callcond)
260         write(*,*) " callcond = ",callcond
261
262         write(*,*)"call thermal conduction in the soil ?"
263         callsoil=.true. ! default value
264         call getin("callsoil",callsoil)
265         write(*,*) " callsoil = ",callsoil
266         
267
268         write(*,*)"call Lott's gravity wave/subgrid topography ",
269     &             "scheme ?"
270         calllott=.true. ! default value
271         call getin("calllott",calllott)
272         write(*,*)" calllott = ",calllott
273
274
275         write(*,*)"rad.transfer is computed every iradia",
276     &             " physical timestep"
277         iradia=1 ! default value
278         call getin("iradia",iradia)
279         write(*,*)" iradia = ",iradia
280         
281
282         write(*,*)"Output of the exchange coefficient mattrix ?",
283     &             "(for diagnostics only)"
284         callg2d=.false. ! default value
285         call getin("callg2d",callg2d)
286         write(*,*)" callg2d = ",callg2d
287
288         write(*,*)"Rayleigh scattering : (should be .false. for now)"
289         rayleigh=.false.
290         call getin("rayleigh",rayleigh)
291         write(*,*)" rayleigh = ",rayleigh
292
293
294! TRACERS:
295
296! dustbin
297         write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)"
298         dustbin=0 ! default value
299         call getin("dustbin",dustbin)
300         write(*,*)" dustbin = ",dustbin
301! active
302         write(*,*)"Radiatively active dust ? (matters if dustbin>0)"
303         active=.false. ! default value
304         call getin("active",active)
305         write(*,*)" active = ",active
306
307! Test of incompatibility:
308! if active is used, then dustbin should be > 0
309
310         if (active.and.(dustbin.lt.1)) then
311           print*,'if active is used, then dustbin should > 0'
312           stop
313         endif
314! doubleq
315         write(*,*)"use mass and number mixing ratios to predict",
316     &             " dust size ?"
317         doubleq=.false. ! default value
318         call getin("doubleq",doubleq)
319         write(*,*)" doubleq = ",doubleq
320! submicron
321         submicron=.false. ! default value
322         call getin("submicron",submicron)
323         write(*,*)" submicron = ",submicron
324
325! Test of incompatibility:
326! if doubleq is used, then dustbin should be 2
327
328         if (doubleq.and.(dustbin.ne.2)) then
329           print*,'if doubleq is used, then dustbin should be 2'
330           stop
331         endif
332         if (doubleq.and.submicron.and.(nq.LT.3)) then
333           print*,'If doubleq is used with a submicron tracer,'
334           print*,' then the number of tracers has to be'
335           print*,' larger than 3.'
336           stop
337         endif
338
339! lifting
340         write(*,*)"dust lifted by GCM surface winds ?"
341         lifting=.false. ! default value
342         call getin("lifting",lifting)
343         write(*,*)" lifting = ",lifting
344
345! Test of incompatibility:
346! if lifting is used, then dustbin should be > 0
347
348         if (lifting.and.(dustbin.lt.1)) then
349           print*,'if lifting is used, then dustbin should > 0'
350           stop
351         endif
352
353! free evolving dust
354! freedust=true just says that there is no lifting and no dust opacity scaling.
355         write(*,*)"dust lifted by GCM surface winds ?"
356         freedust=.false. ! default value
357         call getin("freedust",freedust)
358         write(*,*)" freedust = ",freedust
359         if (freedust.and..not.doubleq) then
360           print*,'freedust should be used with doubleq !'
361           stop
362         endif
363#ifndef MESOSCALE
364         ! this test is valid in GCM case
365         ! ... not in mesoscale case, for we want to activate mesoscale lifting
366         if (freedust.and.lifting) then
367           print*,'if freedust is used, then lifting should not be used'
368           print*,'lifting forced to false !!'
369           lifting=.false.
370         endif
371#endif
372
373! Dust IR opacity
374         write(*,*)" Wavelength for infrared opacity of dust ?"
375         write(*,*)" Choices are:"
376         write(*,*)" tes  --- > 9.3 microns  [default]"
377         write(*,*)" mcs  --- > 21.6 microns"
378         !
379         ! WARNING WARNING WARNING WARNING WARNING WARNING
380         !
381         ! BEFORE ADDING A NEW VALUE, BE SURE THAT THE
382         ! CORRESPONDING WAVELENGTH IS IN THE LOOKUP TABLE,
383         ! OR AT LEAST NO TO FAR, TO AVOID FALLACIOUS INTERPOLATIONS.
384         !
385         dustiropacity="tes" !default value - is expected to shift to mcs one day
386         call getin("dustiropacity",dustiropacity)
387         write(*,*)" dustiropacity = ",trim(dustiropacity)
388         select case (trim(dustiropacity))
389           case ("tes")
390             dustrefir = 9.3E-6
391           case ("mcs")
392             dustrefir = 21.6E-6
393           case default
394              write(*,*) trim(dustiropacity),
395     &                  " is not a valid option for dustiropacity"
396             stop
397         end select
398
399! callddevil
400         write(*,*)" dust lifted by dust devils ?"
401         callddevil=.false. !default value
402         call getin("callddevil",callddevil)
403         write(*,*)" callddevil = ",callddevil
404
405! Test of incompatibility:
406! if dustdevil is used, then dustbin should be > 0
407
408         if (callddevil.and.(dustbin.lt.1)) then
409           print*,'if dustdevil is used, then dustbin should > 0'
410           stop
411         endif
412! sedimentation
413         write(*,*) "Gravitationnal sedimentation ?"
414         sedimentation=.true. ! default value
415         call getin("sedimentation",sedimentation)
416         write(*,*) " sedimentation = ",sedimentation
417! activice
418         write(*,*) "Radiatively active transported atmospheric ",
419     &              "water ice ?"
420         activice=.false. ! default value
421         call getin("activice",activice)
422         write(*,*) " activice = ",activice
423! water
424         write(*,*) "Compute water cycle ?"
425         water=.false. ! default value
426         call getin("water",water)
427         write(*,*) " water = ",water
428
429! thermal inertia feedback
430         write(*,*) "Activate the thermal inertia feedback ?"
431         tifeedback=.false. ! default value
432         call getin("tifeedback",tifeedback)
433         write(*,*) " tifeedback = ",tifeedback
434
435! Test of incompatibility:
436
437         if (tifeedback.and..not.water) then
438           print*,'if tifeedback is used,'
439           print*,'water should be used too'
440           stop
441         endif
442
443         if (tifeedback.and..not.callsoil) then
444           print*,'if tifeedback is used,'
445           print*,'callsoil should be used too'
446           stop
447         endif
448
449         if (activice.and..not.water) then
450           print*,'if activice is used, water should be used too'
451           stop
452         endif
453
454         if (water.and..not.tracer) then
455           print*,'if water is used, tracer should be used too'
456           stop
457         endif
458         
459! water ice clouds effective variance distribution for sedimentaion       
460        write(*,*) "effective variance for water ice clouds ?"
461        nuice_sed=0.45
462        call getin("nuice_sed",nuice_sed)
463        write(*,*) "water_param nueff Sedimentation:", nuice_sed
464         
465! ccn factor if no scavenging         
466        write(*,*) "water param CCN reduc. factor ?", ccn_factor
467        ccn_factor = 4.5
468        call getin("ccn_factor",ccn_factor)
469        write(*,*)" ccn_factor = ",ccn_factor
470        write(*,*)"Careful: only used when microphys=F, otherwise"
471        write(*,*)"the contact parameter is used instead;"
472
473! microphys
474         write(*,*)"Microphysical scheme for water-ice clouds?"
475         microphys=.false. ! default value
476         call getin("microphys",microphys)
477         write(*,*)" microphys = ",microphys
478
479! microphysical parameter contact       
480         write(*,*) "water contact parameter ?"
481         mteta  = 0.95
482         call getin("mteta",mteta)
483         write(*,*) "mteta = ", mteta
484
485! scavenging
486         write(*,*)"Dust scavenging by H2O/CO2 snowfall ?"
487         scavenging=.false. ! default value
488         call getin("scavenging",scavenging)
489         write(*,*)" scavenging = ",scavenging
490         
491
492! Test of incompatibility:
493! if scavenging is used, then dustbin should be > 0
494
495         if ((microphys.and..not.doubleq).or.
496     &       (microphys.and..not.water)) then
497             print*,'if microphys is used, then doubleq,'
498             print*,'and water must be used!'
499             stop
500         endif
501         if (microphys.and..not.scavenging) then
502             print*,''
503             print*,'----------------WARNING-----------------'
504             print*,'microphys is used without scavenging !!!'
505             print*,'----------------WARNING-----------------'
506             print*,''
507         endif
508
509         if ((scavenging.and..not.microphys).or.
510     &       (scavenging.and.(dustbin.lt.1))) then
511             print*,'if scavenging is used, then microphys'
512             print*,'must be used!'
513             stop
514         endif
515
516! Test of incompatibility:
517
518         write(*,*) "Permanent water caps at poles ?",
519     &               " .true. is RECOMMENDED"
520         write(*,*) "(with .true., North cap is a source of water ",
521     &   "and South pole is a cold trap)"
522         caps=.true. ! default value
523         call getin("caps",caps)
524         write(*,*) " caps = ",caps
525
526! albedo_h2o_ice
527         write(*,*) "water ice albedo ?"
528         albedo_h2o_ice=0.45
529         call getin("albedo_h2o_ice",albedo_h2o_ice)
530         write(*,*) " albedo_h2o_ice = ",albedo_h2o_ice
531! inert_h2o_ice
532         write(*,*) "water ice thermal inertia ?"
533         inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
534         call getin("inert_h2o_ice",inert_h2o_ice)
535         write(*,*) " inert_h2o_ice = ",inert_h2o_ice
536! frost_albedo_threshold
537         write(*,*) "frost thickness threshold for albedo ?"
538         frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
539         call getin("frost_albedo_threshold",
540     &    frost_albedo_threshold)
541         write(*,*) " frost_albedo_threshold = ",
542     &            frost_albedo_threshold
543
544! call Titus crocus line -- DEFAULT IS NONE
545         write(*,*) "Titus crocus line ?"
546         tituscap=.false.  ! default value
547         call getin("tituscap",tituscap)
548         write(*,*) "tituscap",tituscap
549                     
550
551         write(*,*) "photochemistry: include chemical species"
552         photochem=.false. ! default value
553         call getin("photochem",photochem)
554         write(*,*) " photochem = ",photochem
555
556
557! SCATTERERS
558         write(*,*) "how many scatterers?"
559         naerkind=1 ! default value
560         call getin("naerkind",naerkind)
561         write(*,*)" naerkind = ",naerkind
562
563! Test of incompatibility
564c        Logical tests for radiatively active water-ice clouds:
565         IF ( (activice.AND.(.NOT.water)).OR.
566     &        (activice.AND.(naerkind.LT.2)) ) THEN
567           WRITE(*,*) 'If activice is TRUE, water has to be set'
568           WRITE(*,*) 'to TRUE, and "naerkind" must be at least'
569           WRITE(*,*) 'equal to 2.'
570           CALL ABORT
571         ELSE IF ( (.NOT.activice).AND.(naerkind.GT.1) ) THEN
572           WRITE(*,*) 'naerkind is greater than unity, but'
573           WRITE(*,*) 'activice has not been set to .true.'
574           WRITE(*,*) 'in callphys.def; this is not logical!'
575           CALL ABORT
576         ENDIF
577
578!------------------------------------------
579!------------------------------------------
580! once naerkind is known allocate arrays
581! -- we do it here and not in phys_var_init
582! -- because we need to know naerkind first
583         CALL ini_scatterers(ngrid,nlayer)
584!------------------------------------------
585!------------------------------------------
586
587
588c        Please name the different scatterers here ----------------
589         name_iaer(1) = "dust_conrath"   !! default choice is good old Conrath profile
590         IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
591         if (nq.gt.1) then
592          ! trick to avoid problems compiling with 1 tracer
593          ! and picky compilers who know name_iaer(2) is out of bounds
594          j=2
595         IF (water.AND.activice) name_iaer(j) = "h2o_ice"      !! radiatively-active clouds
596         IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff
597         endif ! of if (nq.gt.1)
598c        ----------------------------------------------------------
599
600! THERMOSPHERE
601
602         write(*,*) "call thermosphere ?"
603         callthermos=.false. ! default value
604         call getin("callthermos",callthermos)
605         write(*,*) " callthermos = ",callthermos
606         
607
608         write(*,*) " water included without cycle ",
609     &              "(only if water=.false.)"
610         thermoswater=.false. ! default value
611         call getin("thermoswater",thermoswater)
612         write(*,*) " thermoswater = ",thermoswater
613
614         write(*,*) "call thermal conduction ?",
615     &    " (only if callthermos=.true.)"
616         callconduct=.false. ! default value
617         call getin("callconduct",callconduct)
618         write(*,*) " callconduct = ",callconduct
619
620         write(*,*) "call EUV heating ?",
621     &   " (only if callthermos=.true.)"
622         calleuv=.false.  ! default value
623         call getin("calleuv",calleuv)
624         write(*,*) " calleuv = ",calleuv
625
626         write(*,*) "call molecular viscosity ?",
627     &   " (only if callthermos=.true.)"
628         callmolvis=.false. ! default value
629         call getin("callmolvis",callmolvis)
630         write(*,*) " callmolvis = ",callmolvis
631
632         write(*,*) "call molecular diffusion ?",
633     &   " (only if callthermos=.true.)"
634         callmoldiff=.false. ! default value
635         call getin("callmoldiff",callmoldiff)
636         write(*,*) " callmoldiff = ",callmoldiff
637         
638
639         write(*,*) "call thermospheric photochemistry ?",
640     &   " (only if callthermos=.true.)"
641         thermochem=.false. ! default value
642         call getin("thermochem",thermochem)
643         write(*,*) " thermochem = ",thermochem
644
645         write(*,*) "Method to include solar variability"
646         write(*,*) "0-> old method (using solarcondate); ",
647     &                  "1-> variability wit E10.7"
648         solvarmod=1
649         call getin("solvarmod",solvarmod)
650         write(*,*) " solvarmod = ",solvarmod
651
652         write(*,*) "date for solar flux calculation:",
653     &   " (1985 < date < 2002)",
654     $   " (Only used if solvarmod=0)"
655         write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
656         solarcondate=1993.4 ! default value
657         call getin("solarcondate",solarcondate)
658         write(*,*) " solarcondate = ",solarcondate
659         
660         write(*,*) "Solar variability as observed for MY: "
661         write(*,*) "Only if solvarmod=1"
662         solvaryear=24
663         call getin("solvaryear",solvaryear)
664         write(*,*) " solvaryear = ",solvaryear
665
666         write(*,*) "UV heating efficiency:",
667     &   "measured values between 0.19 and 0.23 (Fox et al. 1996)",
668     &   "lower values may be used to compensate low 15 um cooling"
669         euveff=0.21 !default value
670         call getin("euveff",euveff)
671         write(*,*) " euveff = ", euveff
672
673         if (.not.callthermos) then
674           if (thermoswater) then
675             print*,'if thermoswater is set, callthermos must be true'
676             stop
677           endif         
678           if (callconduct) then
679             print*,'if callconduct is set, callthermos must be true'
680             stop
681           endif       
682           if (calleuv) then
683             print*,'if calleuv is set, callthermos must be true'
684             stop
685           endif         
686           if (callmolvis) then
687             print*,'if callmolvis is set, callthermos must be true'
688             stop
689           endif       
690           if (callmoldiff) then
691             print*,'if callmoldiff is set, callthermos must be true'
692             stop
693           endif         
694           if (thermochem) then
695             print*,'if thermochem is set, callthermos must be true'
696             stop
697           endif         
698        endif
699
700! Test of incompatibility:
701! if photochem is used, then water should be used too
702
703         if (photochem.and..not.water) then
704           print*,'if photochem is used, water should be used too'
705           stop
706         endif
707
708! if callthermos is used, then thermoswater should be used too
709! (if water not used already)
710
711         if (callthermos .and. .not.water) then
712           if (callthermos .and. .not.thermoswater) then
713             print*,'if callthermos is used, water or thermoswater
714     &               should be used too'
715             stop
716           endif
717         endif
718
719         PRINT*,'--------------------------------------------'
720         PRINT*
721         PRINT*
722      ELSE
723         write(*,*)
724         write(*,*) 'Cannot read file callphys.def. Is it here ?'
725         stop
726      ENDIF
727
7288000  FORMAT(t5,a12,l8)
7298001  FORMAT(t5,a12,i8)
730
731      PRINT*
732      PRINT*,'conf_phys: daysec',daysec
733      PRINT*
734      PRINT*,'conf_phys: The radiative transfer is computed:'
735      PRINT*,'           each ',iradia,' physical time-step'
736      PRINT*,'        or each ',iradia*dtphys,' seconds'
737      PRINT*
738! --------------------------------------------------------------
739!  Managing the Longwave radiative transfer
740! --------------------------------------------------------------
741
742!     In most cases, the run just use the following values :
743!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
744      callemis=.true.     
745!     ilwd=10*int(daysec/dtphys) ! bug before 22/10/01       
746      ilwd=1
747      ilwn=1 !2
748      ilwb=1 !2
749      linear=.true.       
750      ncouche=3
751      alphan=0.4
752      semi=0
753
754!     BUT people working hard on the LW may want to read them in 'radia.def'
755!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
756
757      OPEN(99,file='radia.def',status='old',form='formatted'
758     .     ,iostat=ierr)
759      IF(ierr.EQ.0) THEN
760         write(*,*) 'conf_phys: Reading radia.def !!!'
761         READ(99,fmt='(a)') ch1
762         READ(99,*) callemis
763         WRITE(*,8000) ch1,callemis
764
765         READ(99,fmt='(a)') ch1
766         READ(99,*) iradia
767         WRITE(*,8001) ch1,iradia
768
769         READ(99,fmt='(a)') ch1
770         READ(99,*) ilwd
771         WRITE(*,8001) ch1,ilwd
772
773         READ(99,fmt='(a)') ch1
774         READ(99,*) ilwn
775         WRITE(*,8001) ch1,ilwn
776
777         READ(99,fmt='(a)') ch1
778         READ(99,*) linear
779         WRITE(*,8000) ch1,linear
780
781         READ(99,fmt='(a)') ch1
782         READ(99,*) ncouche
783         WRITE(*,8001) ch1,ncouche
784
785         READ(99,fmt='(a)') ch1
786         READ(99,*) alphan
787         WRITE(*,*) ch1,alphan
788
789         READ(99,fmt='(a)') ch1
790         READ(99,*) ilwb
791         WRITE(*,8001) ch1,ilwb
792
793
794         READ(99,fmt='(a)') ch1
795         READ(99,'(l1)') callg2d
796         WRITE(*,8000) ch1,callg2d
797
798         READ(99,fmt='(a)') ch1
799         READ(99,*) semi
800         WRITE(*,*) ch1,semi
801      end if
802      CLOSE(99)
803
804      END
Note: See TracBrowser for help on using the repository browser.