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

Last change on this file since 2538 was 2522, checked in by jnaar, 4 years ago

Added flag for temperature-dependant water contact parameter. Default is false for now !
JN

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