Changeset 2304

Timestamp:

Apr 28, 2020, 3:32:34 PM (5 years ago)

Author:

emillour

Message:

Mars GCM:
Some code tidying: use getin_p() instead of getin() and use "call abort_physic"
instead of "stop" or "call abort"
EM

Location:

trunk/LMDZ.MARS

Files:

• README (modified) (1 diff)
• libf/phymars/aeropacity_mod.F (modified) (7 diffs)
• libf/phymars/albedocaps.F90 (modified) (18 diffs)
• libf/phymars/callsedim_mod.F (modified) (8 diffs)
• libf/phymars/conf_phys.F (modified) (82 diffs)
• libf/phymars/improvedclouds_mod.F (modified) (2 diffs)

trunk/LMDZ.MARS/README

@@ -2982 +2982 @@
 Example : concat[diagfi2 (66sols from sol 61); diagfi4 (65sols from sol 193)]
 Old output Time: [61.08,..127,|127.08,..192]; New output Time: [61.08,..127,|193.08,..258]
+
+== 28/04/2020 == EM
+Some code tidying: use getin_p() instead of getin() and use "call abort_physic"
+instead of "stop" or "call abort"

trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F

@@ -12 +12 @@
      &    clearsky,totcloudfrac)
                                                          
-! to use  'getin'
-      USE ioipsl_getincom, only: getin
+      use ioipsl_getin_p_mod, only: getin_p
       use tracer_mod, only: noms, igcm_h2o_ice, igcm_dust_mass,
      &                      igcm_dust_submicron, rho_dust, rho_ice,
@@ -24 +23 @@
       use tracer_mod, only: r3n_q, ref_r0, igcm_dust_number
 #endif
-      use planete_h
-      USE comcstfi_h
+      use comcstfi_h, only: g, pi
       use dimradmars_mod, only: naerkind, name_iaer, 
      &            iaerdust,tauvis,
@@ -32 +30 @@
      &            iaer_stormdust_doubleq,
      &            iaer_topdust_doubleq
-      USE calcstormfract_mod
+
        IMPLICIT NONE
 c=======================================================================
@@ -227 +225 @@
 !       load value of tauvis from callphys.def (if given there,
 !       otherwise default value read from starfi.nc file will be used)
-        call getin("tauvis",tauvis)
+        call getin_p("tauvis",tauvis)
 
         IF (freedust.or.rdstorm) THEN ! if rdstorm no need to held opacity constant at the first levels
@@ -303 +301 @@
         call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
       ELSE
-        stop 'problem with iaervar in aeropacity.F'
+        call abort_physic("aeropacity","wrong value for iaervar",1)
       ENDIF
 
@@ -635 +633 @@
         write(*,*) "Add a local storm ?"
         localstorm=.true. ! default value
-        call getin("localstorm",localstorm)
+        call getin_p("localstorm",localstorm)
         write(*,*) " localstorm = ",localstorm
 
@@ -645 +643 @@
           write(*,*) "ref opacity of local dust storm"
               taulocref = 4.25 ! default value
-              call getin("taulocref",taulocref)
+              call getin_p("taulocref",taulocref)
               write(*,*) " taulocref = ",taulocref
 
           write(*,*) "target altitude of local storm (km)"
               ztoploc = 10.0 ! default value
-              call getin("ztoploc",ztoploc)
+              call getin_p("ztoploc",ztoploc)
               write(*,*) " ztoploc = ",ztoploc
 
           write(*,*) "radius of dust storm (degree)"
               radloc = 0.5 ! default value
-              call getin("radloc",radloc)
+              call getin_p("radloc",radloc)
               write(*,*) " radloc = ",radloc
 
           write(*,*) "center longitude of storm (deg)"
               lonloc = 25.0 ! default value
-              call getin("lonloc",lonloc)
+              call getin_p("lonloc",lonloc)
               write(*,*) " lonloc = ",lonloc
 
           write(*,*) "center latitude of storm (deg)"
               latloc = -2.5 ! default value
-              call getin("latloc",latloc)
+              call getin_p("latloc",latloc)
               write(*,*) " latloc = ",latloc
         
           write(*,*) "reff storm (mic) 0. for background"
               reffstorm = 0.0 ! default value
-              call getin("reffstorm",reffstorm)
+              call getin_p("reffstorm",reffstorm)
               write(*,*) " reffstorm = ",reffstorm
 

trunk/LMDZ.MARS/libf/phymars/albedocaps.F90

@@ -4 +4 @@
 ! depending on the presence of CO2 ice on the surface
 
-! to use the 'getin' routine 
-use ioipsl_getincom, only: getin
+use ioipsl_getin_p_mod, only: getin_p
 use geometry_mod, only: latitude ! grid point latitudes (rad)
 use surfdat_h, only: TESicealbedo, TESice_Ncoef, TESice_Scoef, &
@@ -32 +31 @@
   TESicealbedo=.false. ! default value 
   write(*,*)" albedocaps: Use TES Cap albedoes ?"
-  call getin("TESicealbedo",TESicealbedo)
+  call getin_p("TESicealbedo",TESicealbedo)
   write(*,*)" albedocaps: TESicealbedo = ",TESicealbedo
   
@@ -39 +38 @@
     write(*,*)" albedocaps: Coefficient for Northern Cap ?"
     TESice_Ncoef=1.0 ! default value
-    call getin("TESice_Ncoef",TESice_Ncoef)
+    call getin_p("TESice_Ncoef",TESice_Ncoef)
     write(*,*)" albedocaps: TESice_Ncoef = ",TESice_Ncoef
     
     write(*,*)" albedocaps: Coefficient for Southern Cap ?"
     TESice_Scoef=1.0 ! default value
-    call getin("TESice_Scoef",TESice_Scoef)
+    call getin_p("TESice_Scoef",TESice_Scoef)
     write(*,*)" albedocaps: TESice_Scoef = ",TESice_Scoef
   endif
@@ -116 +115 @@
 real :: lond ! longitude, in degrees
 integer :: i
+character(len=20),parameter :: modname="TES_icecap_albedo"
 
 ! TES datasets: (hard coded fixed length/sizes; for now)
@@ -160 +160 @@
     write(*,*)'   can be obtained online on:'
     write(*,*)'   http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
-    CALL ABORT
+    CALL abort_physic(modname,"missing input file",1)
   ELSE
     write(*,*) "albedocaps: using file ",trim(datadir)//"/npsc_albedo.nc"
@@ -169 +169 @@
     write(*,*) "Failed to find longitude in file!"
     write(*,*)trim(nf90_strerror(ierr))
-    stop
+    call abort_physic(modname,"failed finding longitude",1)
   else
     ierr=nf90_get_var(nid,nvarid,TESlon)
@@ -175 +175 @@
       write(*,*) "Failed loading longitude data from file!"
       write(*,*)trim(nf90_strerror(ierr))
-      stop
+      call abort_physic(modname,"failed loading longitude",1)
     endif
   endif
@@ -183 +183 @@
     write(*,*) "Failed to find latitude in file!"
     write(*,*)trim(nf90_strerror(ierr))
-    stop
+    call abort_physic(modname,"failed finding latitude",1)
   else
     ierr=nf90_get_var(nid,nvarid,TESlatn)
@@ -189 +189 @@
       write(*,*) "Failed loading latitude data from file!"
       write(*,*)trim(nf90_strerror(ierr))
-      stop
+      call abort_physic(modname,"failed loading latitude",1)
     endif
   endif
@@ -197 +197 @@
     write(*,*) "Failed to find time in file!"
     write(*,*)trim(nf90_strerror(ierr))
-    stop
+    call abort_physic(modname,"failed finding time",1)
   else
     ierr=nf90_get_var(nid,nvarid,TESls)
@@ -203 +203 @@
       write(*,*) "Failed loading time data from file!"
       write(*,*)trim(nf90_strerror(ierr))
-      stop
+      call abort_physic(modname,"failed loading time",1)
     endif
   endif
@@ -211 +211 @@
     write(*,*) "Failed to find albedo in file!"
     write(*,*)trim(nf90_strerror(ierr))
-    stop
+    call abort_physic(modname,"failed finding albedo",1)
   else
     ierr=nf90_get_var(nid,nvarid,TESalbn)
@@ -217 +217 @@
       write(*,*) "Failed loading albedo data from file!"
       write(*,*)trim(nf90_strerror(ierr))
-      stop
+      call abort_physic(modname,"failed loading albedo",1)
     endif
   endif
@@ -233 +233 @@
     write(*,*)'   can be obtained online on:'
     write(*,*)'   http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
-    CALL ABORT
+    CALL abort_physic(modname,"missing input file",1)
   ELSE
     write(*,*) "albedocaps: using file ",trim(datadir)//"/spsc_albedo.nc"
@@ -242 +242 @@
     write(*,*) "Failed to find latitude in file!"
     write(*,*)trim(nf90_strerror(ierr))
-    stop
+    call abort_physic(modname,"failed finding latitude",1)
   else
     ierr=nf90_get_var(nid,nvarid,TESlats)
@@ -248 +248 @@
       write(*,*) "Failed loading latitude data from file!"
       write(*,*)trim(nf90_strerror(ierr))
-      stop
+      call abort_physic(modname,"failed loading latitude",1)
     endif
   endif
@@ -256 +256 @@
     write(*,*) "Failed to find albedo in file!"
     write(*,*)trim(nf90_strerror(ierr))
-    stop
+    call abort_physic(modname,"failed finding albedo",1)
   else
     ierr=nf90_get_var(nid,nvarid,TESalbs)
@@ -262 +262 @@
       write(*,*) "Failed loading albedo data from file!"
       write(*,*)trim(nf90_strerror(ierr))
-      stop
+      call abort_physic(modname,"failed loading albedo",1)
     endif
   endif

trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F

@@ -12 +12 @@
      &                tau,tauscaling)
 
-      USE ioipsl_getincom, only: getin
+      USE ioipsl_getin_p_mod, only: getin_p
       USE updaterad, only: updaterdust,updaterice_micro,updaterice_typ
       USE tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
@@ -103 +103 @@
       REAL Mo,No
       REAl ccntyp
-
+      character(len=20),parameter :: modname="callsedim"
 
 
@@ -203 +203 @@
           write(*,*) ' tracers for dust mass and number mixing'
           write(*,*) ' ratio and doubleq is activated!'
-          stop
+          call abort_physic(modname,"missing dust tracers",1)
         endif
         ENDIF !of if (doubleq)
@@ -225 +225 @@
             write(*,*) ' tracers for ccn mass and number mixing'
             write(*,*) ' ratio and microphys is activated!'
-            stop
+            call abort_physic(modname,"missing ccn tracers",1)
           endif
         ENDIF !of if (microphys)
@@ -252 +252 @@
             write(*,*) ' tracers for ccn co2 mass and number mixing'
             write(*,*) ' ratio and co2clouds are activated!'
-            stop
+            call abort_physic(modname,"missing co2 ccn tracers",1)
           endif
        ENDIF                    !of if (co2clouds)
@@ -259 +259 @@
          write(*,*) "correction for the shape of the ice particles ?"
          beta=0.75 ! default value
-         call getin("ice_shape",beta)
+         call getin_p("ice_shape",beta)
          write(*,*) " ice_shape = ",beta
 
@@ -289 +289 @@
              write(*,*) ' tracers for stormdust mass and number mixing'
              write(*,*) ' ratio and rdstorm is activated!'
-             stop
+             call abort_physic(modname,"missing stormdust tracers",1)
            endif
        ENDIF !of if (rdstorm)
@@ -314 +314 @@
              write(*,*) ' tracers for topdust mass and number mixing'
              write(*,*) ' ratio and slpwind is activated!'
-             stop
+             call abort_physic(modname,"missing topdust tracers",1)
            endif
        ENDIF !of if (slpwind)

trunk/LMDZ.MARS/libf/phymars/conf_phys.F

@@ -33 +33 @@
 !   declarations:
 !   -------------
-! to use  'getin'
-      USE ioipsl_getincom, only : getin
       USE ioipsl_getin_p_mod, ONLY : getin_p
       use tracer_mod, only : nuice_sed, ccn_factor, nuiceco2_sed,
@@ -54 +52 @@
       INTEGER,INTENT(IN) :: ngrid,nlayer,nq
       INTEGER ierr,j
+      character(len=20),parameter :: modname="conf_phys"
  
       CHARACTER ch1*12
@@ -60 +59 @@
       ! or shared between dynamics and physics.
       ecritphy=240 ! default value
-      call getin("ecritphy",ecritphy) ! frequency of outputs in physics,
+      call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
                                       ! in dynamical steps
       day_step=960 ! default value
-      call getin("day_step",day_step) ! number of dynamical steps per day
+      call getin_p("day_step",day_step) ! number of dynamical steps per day
       iphysiq=20 ! default value
-      call getin("iphysiq",iphysiq) ! call physics every iphysiq dyn step
+      call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
       ecritstart=0 ! default value
-      call getin("ecritstart",ecritstart) ! write a restart every ecristart steps
+      call getin_p("ecritstart",ecritstart) ! write a restart every ecristart steps
 #endif
 
@@ -88 +87 @@
          write(*,*) "Directory where external input files are:"
          ! default path is set in datafile_mod
-         call getin("datadir",datadir) 
+         call getin_p("datadir",datadir) 
          write(*,*) " datadir = ",trim(datadir)
 
@@ -98 +97 @@
          write(*,*) "Run with or without tracer transport ?"
          tracer=.false. ! default value
-         call getin("tracer",tracer)
+         call getin_p("tracer",tracer)
          write(*,*) " tracer = ",tracer
 
@@ -104 +103 @@
          write(*,*) "(if diurnal=False, diurnal averaged solar heating)"
          diurnal=.true. ! default value
-         call getin("diurnal",diurnal)
+         call getin_p("diurnal",diurnal)
          write(*,*) " diurnal = ",diurnal
 
@@ -111 +110 @@
      &   "set in 'start'"
          season=.true. ! default value
-         call getin("season",season)
+         call getin_p("season",season)
          write(*,*) " season = ",season
 
          write(*,*) "Write some extra output to the screen ?"
          lwrite=.false. ! default value
-         call getin("lwrite",lwrite)
+         call getin_p("lwrite",lwrite)
          write(*,*) " lwrite = ",lwrite
 
@@ -125 +124 @@
          callstats=.true. ! default value
 #endif
-         call getin("callstats",callstats)
+         call getin_p("callstats",callstats)
          write(*,*) " callstats = ",callstats
 
@@ -131 +130 @@
      &              "Climate Database?"
          calleofdump=.false. ! default value
-         call getin("calleofdump",calleofdump)
+         call getin_p("calleofdump",calleofdump)
          write(*,*) " calleofdump = ",calleofdump
 
@@ -139 +138 @@
      &   "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation"
          iaervar=3 ! default value
-         call getin("iaervar",iaervar)
+         call getin_p("iaervar",iaervar)
          write(*,*) " iaervar = ",iaervar
 
@@ -145 +144 @@
      &   "(matters only if iaervar=1)"
          ! NB: default value of tauvis is set/read in startfi.nc file
-         call getin("tauvis",tauvis)
+         call getin_p("tauvis",tauvis)
          write(*,*) " tauvis = ",tauvis
 
@@ -152 +151 @@
      & " =2 Viking scenario; =3 MGS scenario)"
          iddist=3 ! default value
-         call getin("iddist",iddist)
+         call getin_p("iddist",iddist)
          write(*,*) " iddist = ",iddist
 
          write(*,*) "Dust top altitude (km). (Matters only if iddist=1)"
          topdustref= 90.0 ! default value
-         call getin("topdustref",topdustref)
+         call getin_p("topdustref",topdustref)
          write(*,*) " topdustref = ",topdustref
 
          write(*,*) "Prescribed surface thermal flux (H/(rho*cp),K m/s)"
          tke_heat_flux=0. ! default value
-         call getin("tke_heat_flux",tke_heat_flux)
+         call getin_p("tke_heat_flux",tke_heat_flux)
          write(*,*) " tke_heat_flux = ",tke_heat_flux
          write(*,*) " 0 means the usual schemes are computing"
@@ -168 +167 @@
          write(*,*) "call radiative transfer ?"
          callrad=.true. ! default value
-         call getin("callrad",callrad)
+         call getin_p("callrad",callrad)
          write(*,*) " callrad = ",callrad
 
@@ -178 +177 @@
          callslope=.false. ! default value (not supported yet)
 #endif
-         call getin("callslope",callslope)
+         call getin_p("callslope",callslope)
          write(*,*) " callslope = ",callslope
 
@@ -184 +183 @@
      &              "(matters only if callrad=T)"
          callnlte=.false. ! default value
-         call getin("callnlte",callnlte)
+         call getin_p("callnlte",callnlte)
          write(*,*) " callnlte = ",callnlte
          
@@ -193 +192 @@
          write(*,*) "2 -> new model"
          write(*,*) "(matters only if callnlte=T)"
-         call getin("nltemodel",nltemodel)
+         call getin_p("nltemodel",nltemodel)
          write(*,*) " nltemodel = ",nltemodel
 
@@ -199 +198 @@
      &              "(matters only if callrad=T)"
          callnirco2=.false. ! default value
-         call getin("callnirco2",callnirco2)
+         call getin_p("callnirco2",callnirco2)
          write(*,*) " callnirco2 = ",callnirco2
 
@@ -211 +210 @@
          nircorr=0      !default value
 #endif
-         call getin("nircorr",nircorr)
+         call getin_p("nircorr",nircorr)
          write(*,*) " nircorr = ",nircorr
 
          write(*,*) "call turbulent vertical diffusion ?"
          calldifv=.true. ! default value
-         call getin("calldifv",calldifv)
+         call getin_p("calldifv",calldifv)
          write(*,*) " calldifv = ",calldifv
 
          write(*,*) "call thermals ?"
          calltherm=.false. ! default value
-         call getin("calltherm",calltherm)
+         call getin_p("calltherm",calltherm)
          write(*,*) " calltherm = ",calltherm
 
          write(*,*) "call convective adjustment ?"
          calladj=.true. ! default value
-         call getin("calladj",calladj)
+         call getin_p("calladj",calladj)
          write(*,*) " calladj = ",calladj
 
@@ -238 +237 @@
          write(*,*) "used latest version of yamada scheme?"
          callyamada4=.true. ! default value
-         call getin("callyamada4",callyamada4)
+         call getin_p("callyamada4",callyamada4)
          write(*,*) " callyamada4 = ",callyamada4
 
@@ -250 +249 @@
          write(*,*) "call Richardson-based surface layer ?"
          callrichsl=.false. ! default value
-         call getin("callrichsl",callrichsl)
+         call getin_p("callrichsl",callrichsl)
          write(*,*) " callrichsl = ",callrichsl
 
@@ -268 +267 @@
           print*,'If you want to use the Ri. surface-layer, either
      & activate thermals OR de-activate the convective adjustment.'
-          stop
+          call abort_physic(modname,
+     &     "Richardson layer must be used with thermals",1)
          endif
 
          write(*,*) "call CO2 condensation ?"
          callcond=.true. ! default value
-         call getin("callcond",callcond)
+         call getin_p("callcond",callcond)
          write(*,*) " callcond = ",callcond
 
          write(*,*)"call thermal conduction in the soil ?"
          callsoil=.true. ! default value
-         call getin("callsoil",callsoil)
+         call getin_p("callsoil",callsoil)
          write(*,*) " callsoil = ",callsoil
          
@@ -285 +285 @@
      &             "scheme ?"
          calllott=.true. ! default value
-         call getin("calllott",calllott)
+         call getin_p("calllott",calllott)
          write(*,*)" calllott = ",calllott
 
@@ -291 +291 @@
      &             "scheme ?"
          calllott_nonoro=.false. ! default value
-         call getin("calllott_nonoro",calllott_nonoro)
+         call getin_p("calllott_nonoro",calllott_nonoro)
          write(*,*)" calllott_nonoro = ",calllott_nonoro
 
@@ -297 +297 @@
          write(*,*)"call rocket dust storm parametrization"
          rdstorm=.false. ! default value
-         call getin("rdstorm",rdstorm)
+         call getin_p("rdstorm",rdstorm)
          write(*,*)" rdstorm = ",rdstorm
 ! rocket dust storm detrainment coefficient        
         coeff_detrainment=0. ! default value
-        call getin("coeff_detrainment",coeff_detrainment)
+        call getin_p("coeff_detrainment",coeff_detrainment)
         write(*,*)" coeff_detrainment = ",coeff_detrainment
 
@@ -307 +307 @@
          write(*,*)"call slope wind lifting parametrization"
          slpwind=.false. ! default value
-         call getin("slpwind",slpwind)
+         call getin_p("slpwind",slpwind)
          write(*,*)" slpwind = ",slpwind
 
@@ -314 +314 @@
      &              " /condensation of ground water ice"
          latentheat_surfwater=.true. ! default value
-         call getin("latentheat_surfwater",latentheat_surfwater)
+         call getin_p("latentheat_surfwater",latentheat_surfwater)
          write(*,*)" latentheat_surfwater = ",latentheat_surfwater
 
@@ -320 +320 @@
      &             " physical timestep"
          iradia=1 ! default value
-         call getin("iradia",iradia)
+         call getin_p("iradia",iradia)
          write(*,*)" iradia = ",iradia
          
@@ -327 +327 @@
      &             "(for diagnostics only)"
          callg2d=.false. ! default value
-         call getin("callg2d",callg2d)
+         call getin_p("callg2d",callg2d)
          write(*,*)" callg2d = ",callg2d
 
          write(*,*)"Rayleigh scattering : (should be .false. for now)"
          rayleigh=.false.
-         call getin("rayleigh",rayleigh)
+         call getin_p("rayleigh",rayleigh)
          write(*,*)" rayleigh = ",rayleigh
 
@@ -341 +341 @@
          write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)"
          dustbin=0 ! default value
-         call getin("dustbin",dustbin)
+         call getin_p("dustbin",dustbin)
          write(*,*)" dustbin = ",dustbin
 ! active
          write(*,*)"Radiatively active dust ? (matters if dustbin>0)"
          active=.false. ! default value
-         call getin("active",active)
+         call getin_p("active",active)
          write(*,*)" active = ",active
 
@@ -354 +354 @@
          if (active.and.(dustbin.lt.1)) then
            print*,'if active is used, then dustbin should > 0'
-           stop
+           call abort_physic(modname,
+     &          "active option requires dustbin < 0",1)
          endif
 ! doubleq
@@ -360 +361 @@
      &             " dust size ?"
          doubleq=.false. ! default value
-         call getin("doubleq",doubleq)
+         call getin_p("doubleq",doubleq)
          write(*,*)" doubleq = ",doubleq
 ! submicron
          submicron=.false. ! default value
-         call getin("submicron",submicron)
+         call getin_p("submicron",submicron)
          write(*,*)" submicron = ",submicron
 
@@ -372 +373 @@
          if (doubleq.and.(dustbin.ne.2)) then
            print*,'if doubleq is used, then dustbin should be 2'
-           stop
+           call abort_physic(modname,
+     &          "doubleq option requires dustbin = 2",1)
          endif
          if (doubleq.and.submicron.and.(nq.LT.3)) then
@@ -378 +380 @@
            print*,' then the number of tracers has to be'
            print*,' larger than 3.'
-           stop
+           call abort_physic(modname,
+     &          "submicron option requires dustbin > 2",1)
          endif
 
@@ -384 +387 @@
          write(*,*)"dust lifted by GCM surface winds ?"
          lifting=.false. ! default value
-         call getin("lifting",lifting)
+         call getin_p("lifting",lifting)
          write(*,*)" lifting = ",lifting
 
@@ -392 +395 @@
          if (lifting.and.(dustbin.lt.1)) then
            print*,'if lifting is used, then dustbin should > 0'
-           stop
+           call abort_physic(modname,
+     &          "lifting option requires dustbin > 0",1)
          endif
 
 ! dust injection scheme
         dustinjection=0 ! default: no injection scheme
-        call getin("dustinjection",dustinjection)
+        call getin_p("dustinjection",dustinjection)
         write(*,*)" dustinjection = ",dustinjection
 ! dust injection scheme coefficient        
         coeff_injection=1. ! default value
-        call getin("coeff_injection",coeff_injection)
+        call getin_p("coeff_injection",coeff_injection)
         write(*,*)" coeff_in,jection = ",coeff_injection
 ! timing for dust injection        
         ti_injection=10. ! default value
         tf_injection=12. ! default value
-        call getin("ti_injection",ti_injection)
+        call getin_p("ti_injection",ti_injection)
         write(*,*)" ti_injection = ",ti_injection
-        call getin("tf_injection",tf_injection)
+        call getin_p("tf_injection",tf_injection)
         write(*,*)" tf_injection = ",tf_injection
 
@@ -415 +419 @@
          write(*,*)"dust lifted by GCM surface winds ?"
          freedust=.false. ! default value
-         call getin("freedust",freedust)
+         call getin_p("freedust",freedust)
          write(*,*)" freedust = ",freedust
          if (freedust.and..not.doubleq) then
            print*,'freedust should be used with doubleq !'
-           stop
+           call abort_physic(modname,
+     &          "freedust option requires doubleq",1)
          endif
 #ifndef MESOSCALE
@@ -428 +433 @@
              print*,'if freedust is used and dustinjection = 0, 
      &      then lifting should not be used'
-             stop
+             call abort_physic(modname,
+     &          "freedust option with dustinjection = 0"//
+     &          " requires lifting",1)
            endif
          endif
@@ -435 +442 @@
            if(.not.lifting) then
              print*,"if dustinjection=1, then lifting should be true"
-             stop
+             call abort_physic(modname,
+     &          "dustinjection=1 requires lifting",1)
            endif
            if(.not.freedust) then
              print*,"if dustinjection=1, then freedust should be true"
-             stop
+             call abort_physic(modname,
+     &          "dustinjection=1 requires freedust",1)
            endif
          endif
@@ -448 +457 @@
            print*,'if rdstorm or slpwind is used, then doubleq 
      &            should be used !'
-           stop
+           call abort_physic(modname,
+     &          "rdstorm or slpwind requires doubleq",1)
          endif
          if ((rdstorm.or.slpwind).and..not.active) then
            print*,'if rdstorm or slpwind is used, then active 
      &            should be used !'
-           stop
+           call abort_physic(modname,
+     &          "rdstorm or slpwind requires activ",1)
          endif
          if (rdstorm.and..not.lifting) then
            print*,'if rdstorm is used, then lifting 
      &            should be used !'
-           stop
+           call abort_physic(modname,
+     &          "rdstorm requires lifting",1)
          endif
          if ((rdstorm.or.slpwind).and..not.freedust) then
            print*,'if rdstorm or slpwind is used, then freedust 
      &            should be used !'
-           stop
+           call abort_physic(modname,
+     &          "rdstorm or slpwind requires freedust",1)
          endif
          if (rdstorm.and.(dustinjection.eq.0)) then
            print*,'if rdstorm is used, then dustinjection
      &            should be used !'
-           stop
+           call abort_physic(modname,
+     &          "rdstorm requires dustinjection",1)
          endif
 ! Dust IR opacity
@@ -483 +497 @@
          !
          dustiropacity="tes" !default value - is expected to shift to mcs one day
-         call getin("dustiropacity",dustiropacity)
+         call getin_p("dustiropacity",dustiropacity)
          write(*,*)" dustiropacity = ",trim(dustiropacity)
          select case (trim(dustiropacity))
@@ -493 +507 @@
               write(*,*) trim(dustiropacity),
      &                  " is not a valid option for dustiropacity"
-             stop
+             call abort_physic(modname,
+     &          "invalid dustiropacity option value",1)
          end select
 
@@ -499 +514 @@
          write(*,*)" dust lifted by dust devils ?"
          callddevil=.false. !default value
-         call getin("callddevil",callddevil)
+         call getin_p("callddevil",callddevil)
          write(*,*)" callddevil = ",callddevil
 
@@ -507 +522 @@
          if (callddevil.and.(dustbin.lt.1)) then
            print*,'if dustdevil is used, then dustbin should > 0'
-           stop
+           call abort_physic(modname,
+     &          "callddevil requires dustbin > 0",1)
          endif
 ! sedimentation
          write(*,*) "Gravitationnal sedimentation ?"
          sedimentation=.true. ! default value
-         call getin("sedimentation",sedimentation)
+         call getin_p("sedimentation",sedimentation)
          write(*,*) " sedimentation = ",sedimentation
 ! activice
@@ -518 +534 @@
      &              "water ice ?"
          activice=.false. ! default value
-         call getin("activice",activice)
+         call getin_p("activice",activice)
          write(*,*) " activice = ",activice
 ! water
          write(*,*) "Compute water cycle ?"
          water=.false. ! default value
-         call getin("water",water)
+         call getin_p("water",water)
          write(*,*) " water = ",water
 
@@ -529 +545 @@
          write(*,*) "Fixed cloud fraction?"
          CLFfixval=1.0 ! default value
-         call getin("CLFfixval",CLFfixval)
+         call getin_p("CLFfixval",CLFfixval)
          write(*,*) "CLFfixval=",CLFfixval
 ! sub-grid cloud fraction: varying
          write(*,*) "Use partial nebulosity?"
          CLFvarying=.false. ! default value
-         call getin("CLFvarying",CLFvarying)
+         call getin_p("CLFvarying",CLFvarying)
          write(*,*)"CLFvarying=",CLFvarying
 
@@ -540 +556 @@
          write(*,*) "Compute CO2 clouds (implies microphysical scheme)?"
          co2clouds=.false. ! default value
-         call getin("co2clouds",co2clouds)
+         call getin_p("co2clouds",co2clouds)
          write(*,*) " co2clouds = ",co2clouds
 !Can water ice particles serve as CCN for CO2clouds
          write(*,*) "Use water ice as CO2 clouds CCN ?"
          co2useh2o=.false. ! default value
-         call getin("co2useh2o",co2useh2o)
+         call getin_p("co2useh2o",co2useh2o)
          write(*,*) " co2useh2o = ",co2useh2o
 !Do we allow a supply of meteoritic paricles to serve as CO2 ice CCN?
          write(*,*) "Supply meteoritic particle for CO2 clouds ?"
          meteo_flux=.false. !Default value
-         call getin("meteo_flux",meteo_flux)
+         call getin_p("meteo_flux",meteo_flux)
          write(*,*)  " meteo_flux = ",meteo_flux
 !Do we allow a sub-grid temperature distribution for the CO2 microphysics
          write(*,*) "sub-grid temperature distribution for CO2 clouds?"
          CLFvaryingCO2=.false. !Default value
-         call getin("CLFvaryingCO2",CLFvaryingCO2)
+         call getin_p("CLFvaryingCO2",CLFvaryingCO2)
          write(*,*)  " CLFvaryingCO2 = ",CLFvaryingCO2
 !Amplitude of the sub-grid temperature distribution for the CO2 microphysics
          write(*,*) "sub-grid temperature amplitude for CO2 clouds?"
          spantCO2=0 !Default value
-         call getin("spantCO2",spantCO2)
+         call getin_p("spantCO2",spantCO2)
          write(*,*)  " spantCO2 = ",spantCO2
 !Do you want to filter the sub-grid T distribution by a Saturation index?
          write(*,*) "filter sub-grid temperature by Saturation index?"
          satindexco2=.true.
-         call getin("satindexco2",satindexco2)
+         call getin_p("satindexco2",satindexco2)
          write(*,*)  " satindexco2 = ",satindexco2
 
@@ -572 +588 @@
          write(*,*) "Activate the thermal inertia feedback ?"
          tifeedback=.false. ! default value
-         call getin("tifeedback",tifeedback)
+         call getin_p("tifeedback",tifeedback)
          write(*,*) " tifeedback = ",tifeedback
 
@@ -580 +596 @@
            print*,'if tifeedback is used,'
            print*,'water should be used too'
-           stop
+           call abort_physic(modname,
+     &          "tifeedback requires water",1)
          endif
 
@@ -586 +603 @@
            print*,'if tifeedback is used,'
            print*,'callsoil should be used too'
-           stop
+           call abort_physic(modname,
+     &          "tifeedback requires callsoil",1)
          endif
 
          if (activice.and..not.water) then
            print*,'if activice is used, water should be used too'
-           stop
+           call abort_physic(modname,
+     &          "activeice requires water",1)
          endif
 
          if (water.and..not.tracer) then
            print*,'if water is used, tracer should be used too'
-           stop
+           call abort_physic(modname,
+     &          "water requires tracer",1)
          endif
          
@@ -602 +622 @@
         write(*,*) "Sed effective variance for water ice clouds ?"
         nuice_sed=0.45 
-        call getin("nuice_sed",nuice_sed)
+        call getin_p("nuice_sed",nuice_sed)
         write(*,*) "water_param nueff Sedimentation:", nuice_sed
               
         write(*,*) "Sed effective variance for CO2 clouds ?"
         nuiceco2_sed=0.45 
-        call getin("nuiceco2_sed",nuiceco2_sed)
+        call getin_p("nuiceco2_sed",nuiceco2_sed)
         write(*,*) "CO2 nueff Sedimentation:", nuiceco2_sed
   
         write(*,*) "REF effective variance for CO2 clouds ?"
         nuiceco2_ref=0.45 
-        call getin("nuiceco2_ref",nuiceco2_ref)
+        call getin_p("nuiceco2_ref",nuiceco2_ref)
         write(*,*) "CO2 nueff Sedimentation:", nuiceco2_ref
 
         write(*,*) "REF effective variance for water clouds ?"
         nuice_ref=0.45 
-        call getin("nuice_ref",nuice_ref)
+        call getin_p("nuice_ref",nuice_ref)
         write(*,*) "CO2 nueff Sedimentation:", nuice_ref
 
@@ -624 +644 @@
         write(*,*) "water param CCN reduc. factor ?"
         ccn_factor = 4.5
-        call getin("ccn_factor",ccn_factor)
+        call getin_p("ccn_factor",ccn_factor)
         write(*,*)" ccn_factor = ",ccn_factor
         write(*,*)"Careful: only used when microphys=F, otherwise"
@@ -632 +652 @@
         write(*,*)"Microphysical scheme for water-ice clouds?"
         microphys=.false.       ! default value
-        call getin("microphys",microphys)
+        call getin_p("microphys",microphys)
         write(*,*)" microphys = ",microphys
 
@@ -638 +658 @@
         write(*,*)"Allow super-saturation of water vapor?"
         supersat=.true.         ! default value
-        call getin("supersat",supersat)
+        call getin_p("supersat",supersat)
         write(*,*)"supersat = ",supersat
 
@@ -644 +664 @@
         write(*,*) "water contact parameter ?"
         mteta  = 0.95
-        call getin("mteta",mteta)
+        call getin_p("mteta",mteta)
         write(*,*) "mteta = ", mteta
         
@@ -650 +670 @@
         write(*,*)"Dust scavenging by H2O/CO2 snowfall ?"
         scavenging=.false.      ! default value
-        call getin("scavenging",scavenging)
+        call getin_p("scavenging",scavenging)
         write(*,*)" scavenging = ",scavenging
          
@@ -661 +681 @@
            print*,'if microphys is used, then doubleq,'
            print*,'and water must be used!'
-           stop
+           call abort_physic(modname,
+     &          "microphys requires water and doubleq",1)
         endif
         if (microphys.and..not.scavenging) then
@@ -675 +696 @@
            print*,'if scavenging is used, then microphys'
            print*,'must be used!'
-           stop
+           call abort_physic(modname,
+     &          "scavenging requires microphys",1)
         endif
 
@@ -682 +704 @@
         write(*,*)"Dust scavenging by instantaneous CO2 snowfall ?"
         scavco2cond=.false.      ! default value
-        call getin("scavco2cond",scavco2cond)
+        call getin_p("scavco2cond",scavco2cond)
         write(*,*)" scavco2cond = ",scavco2cond
 ! Test of incompatibility:
@@ -688 +710 @@
         if (scavco2cond.and.(dustbin.lt.1))then
            print*,'if scavco2cond is used, then dustbin should be > 0'
-           stop
+           call abort_physic(modname,
+     &          "scavco2cond requires dustbin > 0",1)
         endif
 ! if co2clouds is used, then there is no need for scavco2cond
@@ -698 +721 @@
            print*,'----------------WARNING-----------------'
            print*,''
-           stop
+           call abort_physic(modname,
+     &          "incompatible co2cloud and scavco2cond options",1)
         endif
         
@@ -708 +732 @@
      &   "and South pole is a cold trap)"
          caps=.true. ! default value
-         call getin("caps",caps)
+         call getin_p("caps",caps)
          write(*,*) " caps = ",caps
 
@@ -714 +738 @@
          write(*,*) "water ice albedo ?"
          albedo_h2o_ice=0.45
-         call getin("albedo_h2o_ice",albedo_h2o_ice)
+         call getin_p("albedo_h2o_ice",albedo_h2o_ice)
          write(*,*) " albedo_h2o_ice = ",albedo_h2o_ice
 ! inert_h2o_ice
          write(*,*) "water ice thermal inertia ?"
          inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
-         call getin("inert_h2o_ice",inert_h2o_ice)
+         call getin_p("inert_h2o_ice",inert_h2o_ice)
          write(*,*) " inert_h2o_ice = ",inert_h2o_ice
 ! frost_albedo_threshold
          write(*,*) "frost thickness threshold for albedo ?"
          frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
-         call getin("frost_albedo_threshold",
+         call getin_p("frost_albedo_threshold",
      &    frost_albedo_threshold)
          write(*,*) " frost_albedo_threshold = ",
@@ -732 +756 @@
          write(*,*) "Titus crocus line ?"
          tituscap=.false.  ! default value
-         call getin("tituscap",tituscap)
+         call getin_p("tituscap",tituscap)
          write(*,*) "tituscap",tituscap
                      
@@ -738 +762 @@
          write(*,*) "photochemistry: include chemical species"
          photochem=.false. ! default value
-         call getin("photochem",photochem)
+         call getin_p("photochem",photochem)
          write(*,*) " photochem = ",photochem
          
@@ -744 +768 @@
      &   "every ichemistry physics step (default: ichemistry=1)"
          ichemistry=1
-         call getin("ichemistry",ichemistry)
+         call getin_p("ichemistry",ichemistry)
          write(*,*) " ichemistry = ",ichemistry
 
@@ -751 +775 @@
          write(*,*) "how many scatterers?"
          naerkind=1 ! default value
-         call getin("naerkind",naerkind)
+         call getin_p("naerkind",naerkind)
          write(*,*)" naerkind = ",naerkind
 
@@ -761 +785 @@
            WRITE(*,*) 'to TRUE, and "naerkind" must be at least'
            WRITE(*,*) 'equal to 2.'
-           CALL ABORT
+           call abort_physic(modname,
+     &          "radiatively active dust and water"//
+     &          " require naerkind > 1",1)
          ENDIF
 
@@ -807 +833 @@
          write(*,*) "call thermosphere ?"
          callthermos=.false. ! default value
-         call getin("callthermos",callthermos)
+         call getin_p("callthermos",callthermos)
          write(*,*) " callthermos = ",callthermos
          
@@ -814 +840 @@
      &              "(only if water=.false.)"
          thermoswater=.false. ! default value
-         call getin("thermoswater",thermoswater)
+         call getin_p("thermoswater",thermoswater)
          write(*,*) " thermoswater = ",thermoswater
 
@@ -820 +846 @@
      &    " (only if callthermos=.true.)"
          callconduct=.false. ! default value
-         call getin("callconduct",callconduct)
+         call getin_p("callconduct",callconduct)
          write(*,*) " callconduct = ",callconduct
 
@@ -826 +852 @@
      &   " (only if callthermos=.true.)"
          calleuv=.false.  ! default value
-         call getin("calleuv",calleuv)
+         call getin_p("calleuv",calleuv)
          write(*,*) " calleuv = ",calleuv
 
@@ -832 +858 @@
      &   " (only if callthermos=.true.)"
          callmolvis=.false. ! default value
-         call getin("callmolvis",callmolvis)
+         call getin_p("callmolvis",callmolvis)
          write(*,*) " callmolvis = ",callmolvis
 
@@ -838 +864 @@
      &   " (only if callthermos=.true.)"
          callmoldiff=.false. ! default value
-         call getin("callmoldiff",callmoldiff)
+         call getin_p("callmoldiff",callmoldiff)
          write(*,*) " callmoldiff = ",callmoldiff
          
@@ -845 +871 @@
      &   " (only if callthermos=.true.)"
          thermochem=.false. ! default value
-         call getin("thermochem",thermochem)
+         call getin_p("thermochem",thermochem)
          write(*,*) " thermochem = ",thermochem
 
@@ -852 +878 @@
      &          "1-> daily evolution of E10.7 (for given solvaryear)"
          solvarmod=1
-         call getin("solvarmod",solvarmod)
+         call getin_p("solvarmod",solvarmod)
          write(*,*) " solvarmod = ",solvarmod
 
@@ -858 +884 @@
          write(*,*) " (min=80 , ave=140, max=320)"
          fixed_euv_value=140 ! default value
-         call getin("fixed_euv_value",fixed_euv_value)
+         call getin_p("fixed_euv_value",fixed_euv_value)
          write(*,*) " fixed_euv_value = ",fixed_euv_value
          
@@ -864 +890 @@
          write(*,*) "Only if solvarmod=1"
          solvaryear=24
-         call getin("solvaryear",solvaryear)
+         call getin_p("solvaryear",solvaryear)
          write(*,*) " solvaryear = ",solvaryear
 
@@ -871 +897 @@
      &   "lower values may be used to compensate low 15 um cooling"
          euveff=0.21 !default value
-         call getin("euveff",euveff)
+         call getin_p("euveff",euveff)
          write(*,*) " euveff = ", euveff
 
@@ -878 +904 @@
            if (thermoswater) then
              print*,'if thermoswater is set, callthermos must be true'
-             stop
+             call abort_physic(modname,
+     &          "thermoswater requires callthermos",1)
            endif          
            if (callconduct) then
              print*,'if callconduct is set, callthermos must be true'
-             stop
+             call abort_physic(modname,
+     &          "callconduct requires callthermos",1)
            endif        
            if (calleuv) then
              print*,'if calleuv is set, callthermos must be true'
-             stop
+             call abort_physic(modname,
+     &          "calleuv requires callthermos",1)
            endif         
            if (callmolvis) then
              print*,'if callmolvis is set, callthermos must be true'
-             stop
+             call abort_physic(modname,
+     &          "callmolvis requires callthermos",1)
            endif        
            if (callmoldiff) then
              print*,'if callmoldiff is set, callthermos must be true'
-             stop
+             call abort_physic(modname,
+     &          "callmoldiff requires callthermos",1)
            endif          
            if (thermochem) then
              print*,'if thermochem is set, callthermos must be true'
-             stop
+             call abort_physic(modname,
+     &          "thermochem requires callthermos",1)
            endif          
         endif
@@ -907 +939 @@
          if (photochem.and..not.water) then
            print*,'if photochem is used, water should be used too'
-           stop
+           call abort_physic(modname,
+     &          "photochem requires water",1)
          endif
 
@@ -917 +950 @@
              print*,'if callthermos is used, water or thermoswater 
      &               should be used too'
-             stop
+             call abort_physic(modname,
+     &          "callthermos requires water or thermoswater",1)
            endif
          endif
@@ -927 +961 @@
          write(*,*)
          write(*,*) 'Cannot read file callphys.def. Is it here ?'
-         stop
+         call abort_physic(modname,
+     &          "missing callphys.def file",1)
       ENDIF
 

trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F

@@ -9 +9 @@
      &             pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
      &             nq,tauscaling)
-! to use  'getin'
-      USE ioipsl_getincom
-      USE updaterad
+      USE updaterad, ONLY: updaterice_micro, updaterccn
       USE watersat_mod, ONLY: watersat
       use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap,
@@ -18 +16 @@
      &                      igcm_ccn_number
       use conc_mod, only: mmean
-      USE comcstfi_h
+      use comcstfi_h, only: pi, cpp
       implicit none