Changeset 120

Timestamp:

May 17, 2011, 4:27:01 PM (14 years ago)

Author:

emillour

Message:

-Mars GCM:

set internal computations using double precision in growthrate.F and

watercloud.F (otherwise we sometimes end up with Nans).

add extra checks in newcondens.F to avoid possibility of out of bounds

evaluation of array masse()

EM

Location:

trunk/mars

Files:

• README (modified) (1 diff)
• libf/phymars/growthrate.F (modified) (3 diffs)
• libf/phymars/newcondens.F (modified) (1 diff)
• libf/phymars/newsedim.F (modified) (1 diff)
• libf/phymars/watercloud.F (modified) (5 diffs)

trunk/mars/README

@@ -632 +632 @@
    --> fix by moving the closing "endif" [in addition to corrections mentioned in the previous point]
 
+== 17/10/2011 == EM
+>> set internal computations using double precision in growthrate.F and
+   watercloud.F (otherwise we sometimes end up with Nans).
+>> add extra checks in newcondens.F to avoid possibility of out of bounds
+   evaluation of array masse()
+   

trunk/mars/libf/phymars/growthrate.F

@@ -23 +23 @@
       REAL*8 psat ! water vapor saturation pressure (Pa) 
       REAL r    ! crystal radius before condensation (m)
-      REAL seq  ! Equilibrium saturation ratio
-      REAL dr   ! crystal radius variation (m)
+      REAL*8 seq  ! Equilibrium saturation ratio
+!      REAL dr   ! crystal radius variation (m)
+      REAL*8 Cste 
 
 c   local:
 c   ------
 
-      REAL molco2,molh2o
-      REAL Mco2,Mh2o,rho_i,sigh2o
-      REAL nav,rgp,kbz,pi,To
+      REAL*8 molco2,molh2o
+      REAL*8 Mco2,Mh2o,rho_i,sigh2o
+      REAL*8 nav,rgp,kbz,pi,To
 
 c     Effective gas molecular radius (m)
-      data molco2/2.2e-10/   ! CO2
+      data molco2/2.2d-10/   ! CO2
 c     Effective gas molecular radius (m)
-      data molh2o/1.2e-10/   ! H2O
+      data molh2o/1.2d-10/   ! H2O
 c     Molecular weight of CO2
-      data Mco2/44.e-3/            ! kg.mol-1
+      data Mco2/44.d-3/            ! kg.mol-1
 c     Molecular weight of H2O
-      data Mh2o/18.e-3/            ! kg.mol-1
+      data Mh2o/18.d-3/            ! kg.mol-1
 c     surface tension of ice/vapor
       data sigh2o/0.12/      ! N.m
@@ -46 +47 @@
       data rho_i/917./        ! kg.m-3 also defined in initcld.f
 c     Avogadro number
-      data nav/6.023e23/
+      data nav/6.023d23/
 c     Perfect gas constant
       data rgp/8.3143/
 c     Boltzman constant
-      data kbz/1.381e-23/
+      data kbz/1.381d-23/
 c     pi number
       data pi/3.141592654/
@@ -56 +57 @@
       data To/273.15/
 
-      REAL k,Lv                 
-      REAL knudsen           ! Knudsen number (gas mean free path/particle radius)
-      REAL a,Dv,lambda       ! Intermediate computations for growth rate
+      REAL*8 k,Lv                 
+      REAL*8 knudsen           ! Knudsen number (gas mean free path/particle radius)
+      REAL*8 a,Dv,lambda       ! Intermediate computations for growth rate
       REAL*8 Rk,Rd
-      REAL Cste, rf
 
 c-----------------------------------------------------------------------

trunk/mars/libf/phymars/newcondens.F

@@ -854 +854 @@
              Mtot = masse(m+1)
              MQtot = masse(m+1)*q(m+1)
-             do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
+             if (m.gt.0) then ! because some compilers will have problems
+                              ! evaluating masse(0)
+              do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
                 m=m-1
                 Mtot = Mtot + masse(m+1)
                 MQtot = MQtot + masse(m+1)*q(m+1)
-             end do
+                if (m.eq.0) exit
+              end do
+             endif
              if (m.gt.0) then
                 sigw=(w(l+1)+Mtot)/masse(m)

trunk/mars/libf/phymars/newsedim.F

@@ -179 +179 @@
              ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
              IF ( ptop .eq. 1. ) THEN
-                PRINT*, 'exposant trop petit ', ig, l
+                PRINT*, 'newsedim: exposant trop petit ', ig, l
                 ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k)))
              ELSE

trunk/mars/libf/phymars/watercloud.F

@@ -75 +75 @@
       REAL CBRT
       EXTERNAL CBRT
-      INTEGER ig,iq,l
+      INTEGER ig,l
 
 
@@ -85 +85 @@
       REAL masse (ngridmx,nlayermx)
       REAL epaisseur (ngridmx,nlayermx)
-      REAL rfinal        ! Ice crystal radius after condensation(m)
-      REAL seq           ! Equilibrium saturation ration (accounting for curvature effect)
-      REAL dzq           ! masse de glace echangee (kg/kg)
+!      REAL rfinal        ! Ice crystal radius after condensation(m)
+      REAL*8 seq           ! Equilibrium saturation ration (accounting for curvature effect)
+      REAL*8 dzq           ! masse de glace echangee (kg/kg)
       REAL lw       !Latent heat of sublimation (J.kg-1) 
       REAL,PARAMETER :: To=273.15 ! reference temperature, T=273.15 K
@@ -93 +93 @@
       REAL Ctot
       REAL*8 ph2o,satu
-      REAL gr,Cste,up,dwn,newvap
+      REAL*8 gr,Cste,up,dwn,newvap
 
       LOGICAL,SAVE :: firstcall=.true.
@@ -104 +104 @@
 c     REAL rave(ngridmx)               ! Mean crystal radius in a column (m)
 
-      INTEGER i
-
 ! indexes of water vapour, water ice and dust tracers:
       INTEGER,SAVE :: i_h2o=0  ! water vapour
       INTEGER,SAVE :: i_ice=0  ! water ice
-      CHARACTER(LEN=20) :: tracername ! to temporarly store text
 
 c    ** un petit test de coherence
@@ -209 +206 @@
 c         Improved microphysics scheme
 c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
            Ctot = zq(ig,l,i_h2o) + zq(ig,l,i_ice)
            ph2o = zq(ig,l,i_h2o) * 44. / 18. * pplay(ig,l)