Changeset 1884


Ignore:
Timestamp:
Jan 5, 2018, 1:27:35 PM (7 years ago)
Author:
emillour
Message:

Mars GCM:

  • updated massflowrateCO2 routine which now uses an analytical formula rather than an iterative solver
  • some code tidying in improvedCO2clouds.F

CL+EM

Location:
trunk/LMDZ.MARS
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r1881 r1884  
    24952495-zteff and pqeff are initialized in the first part of the CLFvarying scheme, section 0 of the code, instead of been initialized in section 1 (tendencies) with the sub-timesteps.
    24962496-The cloud fraction cannot be lower than the settled value mincloud.
     2497
     2498== 05/01/2018 == CL+EM
     2499- updated massflowrateCO2 routine which now uses an analytical formula rather
     2500  than an iterative solver
     2501- some code tidying in improvedCO2clouds.F

  • trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F

    r1818 r1884  
    44     &             nq,tauscaling,
    55     &             memdMMccn,memdMMh2o,memdNNccn)
    6       USE comcstfi_h
    7       USE ioipsl_getincom
    8       USE updaterad
    9       use tracer_mod
     6      USE comcstfi_h, only: pi, g, cpp
     7      USE updaterad, only: updaterice_micro, updaterice_microco2
     8      use tracer_mod, only: igcm_dust_mass, igcm_dust_number, rho_dust,
     9     &                      igcm_h2o_ice, igcm_ccn_mass,
     10     &                      igcm_ccn_number, nuice_sed,
     11     &                      igcm_co2, igcm_co2_ice, igcm_ccnco2_mass,
     12     &                      igcm_ccnco2_number, nuiceco2_sed,
     13     &                      rho_ice_co2
    1014      use conc_mod, only: mmean
    1115
    … …  
    4953 
    5054c------------------------------------------------------------------
    51 #include "callkeys.h"
    52 #include "microphys.h"
    53 #include "datafile.h"
     55      include "callkeys.h"
     56      include "microphys.h"
     57      include "datafile.h"
    5458c------------------------------------------------------------------
    55 c     Inputs:
    56 
    57       INTEGER ngrid,nlay
    58       integer nq                 ! nombre de traceurs
    59       REAL ptimestep             ! pas de temps physique (s)
    60       REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
    61       REAL pplev(ngrid,nlay+1)   ! pression inter couches (Pa)
     59c     Arguments:
     60
     61      INTEGER,INTENT(in) :: ngrid,nlay
     62      integer,intent(in) :: nq         ! number of tracers
     63      REAL,INTENT(in) :: ptimestep     ! physics time step (s)
     64      REAL,INTENT(in) :: pplay(ngrid,nlay)     ! mid-layer pressure (Pa)
     65      REAL,INTENT(in) :: pplev(ngrid,nlay+1)   ! inter-layer pressure (Pa)
     66      REAL,INTENT(in) :: pt(ngrid,nlay) ! temperature at the middle of the
     67                                 !   layers (K)
     68      REAL,INTENT(in) :: pdt(ngrid,nlay) ! tendency on temperature from
     69                                 !  previous physical parametrizations
     70      REAL,INTENT(in) :: pq(ngrid,nlay,nq) ! tracers (kg/kg)
     71      REAL,INTENT(in) :: pdq(ngrid,nlay,nq) ! tendencies on tracers
     72                                 !  before condensation (kg/kg.s-1)
     73      REAL,INTENT(in) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust
     74c     Outputs:
     75      REAL,INTENT(out) :: pdqcloudco2(ngrid,nlay,nq) ! tendency on tracers
     76                                   ! due to CO2 condensation (kg/kg.s-1)
     77      ! condensation si igcm_co2_ice
     78      REAL,INTENT(out) :: pdtcloudco2(ngrid,nlay)  ! tendency on temperature due
     79                                   ! to latent heat
     80
     81c------------------------------------------------------------------
     82c     Local variables:
     83      LOGICAL,SAVE :: firstcall=.true.
     84      REAL*8   derf ! Error function
     85      INTEGER ig,l,i,flag_pourri
     86     
    6287      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
    63       REAL pt(ngrid,nlay)        ! temperature at the middle of the
    64                                  !   layers (K)
    65       REAL pdt(ngrid,nlay)       ! tendance temperature des autres
    66                                  !   param.
    67       REAL pq(ngrid,nlay,nq)     ! traceur (kg/kg)
    68       REAL pdq(ngrid,nlay,nq)    ! tendance avant condensation
    69                                  !   (kg/kg.s-1)
    70       REAL tauscaling(ngrid)     ! Convertion factor for qdust and Ndust
    7188      REAL rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
    7289                                ! used for nucleation of CO2 on ice-coated ccns
    7390      REAL rccnh2o(ngrid,nlay)    ! Water Ice mass mean radius (m)
    74 c     Outputs:
    75       REAL pdqcloudco2(ngrid,nlay,nq) ! tendance de la condensation
    76                                    !   CO2 (kg/kg.s-1)
    77       ! condensation si igcm_co2_ice
    78       REAL pdtcloudco2(ngrid,nlay)    ! tendance temperature due
    79                                    !   a la chaleur latente
    80 
    81 c------------------------------------------------------------------
    82 c     Local variables:
    83       LOGICAL firstcall
    84       DATA firstcall/.true./
    85       SAVE firstcall
    86       REAL*8   derf ! Error function
    87       INTEGER ig,l,i,flag_pourri
    88      
    8991      REAL zq(ngrid,nlay,nq)  ! local value of tracers
    9092      REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
    … …  
    156158      integer,save :: coeffh2o  !will be =0 is co2useh2o=.false.
    157159!     Variables for the meteoritic flux:
    158       double precision meteo_ccn(ngrid,nlay,100) !100=nbinco2_cld !!!
    159       double precision,save :: meteo(130,100)
    160       double precision mtemp(100),pression_meteo(130)
     160      integer,parameter :: nbin_meteor=100
     161      integer,parameter :: nlev_meteor=130
     162      double precision meteor_ccn(ngrid,nlay,100) !100=nbinco2_cld !!!
     163      double precision,save :: meteor(130,100)
     164      double precision mtemp(100),pression_meteor(130)
    161165      logical file_ok
     166      integer read_ok
    162167      integer nelem,lebon1,lebon2
    163168      double precision :: ltemp1(130),ltemp2(130)
    164       integer ibin,uMeteo,j
     169      integer ibin,uMeteor,j
    165170
    166171      IF (firstcall) THEN
    … …  
    198203        print*,'-----------------------------------'
    199204        do i=1,nbinco2_cld
    200           write(*,'(i3,3x,3(e12.6,4x))') i,rb_cldco2(i), rad_cldco2(i),
     205          write(*,'(i3,3x,3(e13.6,4x))') i,rb_cldco2(i), rad_cldco2(i),
    201206     &      dr_cld(i)
    202207        enddo
    203         write(*,'(i3,3x,e12.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1)
     208        write(*,'(i3,3x,e13.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1)
    204209        print*,'-----------------------------------'
    205210        do i=1,nbinco2_cld+1
    … …  
    249254           endif
    250255!used Variables
    251            open(newunit=uMeteo,file=trim(datafile)//
     256           open(newunit=uMeteor,file=trim(datafile)//
    252257     &          '/Meteo_flux_Plane.dat'
    253258     &          ,FORM='formatted')
    254259!13000 records (130 pressions x 100 bin sizes)
    255            read(uMeteo,*) !skip 1 line
     260           read(uMeteor,*) !skip 1 line
    256261           do i=1,130
    257               read(uMeteo,*) pression_meteo(i)
    258               write(*,*) pression_meteo(i)
     262              read(uMeteor,*,iostat=read_ok) pression_meteor(i)
     263              if (read_ok==0) then
     264                write(*,*) pression_meteor(i)
     265              else
     266                write(*,*) 'Error reading Meteo_flux_Plane.dat'
     267                call abort_physic("CO2clouds",
     268     &                    "Error reading Meteo_flux_Plane.dat",1)
     269              endif
    259270           enddo
    260            read(uMeteo,*)       !skip 1 line
     271           read(uMeteor,*)       !skip 1 line
    261272           do i=1,130
    262273              do j=1,100        ! les mêmes 100 bins size que la distri nuclea: on touche pas
    263                  read(uMeteo,'(F12.6)') meteo(i,j)
     274                 read(uMeteor,'(F12.6)',iostat=read_ok) meteor(i,j)
     275                 if (read_ok/=0) then
     276                   write(*,*) 'Error reading Meteo_flux_Plane.dat'
     277                   call abort_physic("CO2clouds",
     278     &                    "Error reading Meteo_flux_Plane.dat",1)
     279                 endif
    264280              enddo
    265281!On doit maintenant réinterpoler le tableau(130,100) sur les pressions du GCM (nlay,100)
    266282           enddo
    267            close(uMeteo)
     283           close(uMeteor)
    268284        write(*,*) "File meteo_flux read, end of firstcall in co2 micro"
    269285        endif                     !of if meteo_flux
    … …  
    281297!=============================================================
    282298      flag_pourri=0
    283       meteo_ccn(:,:,:)=0.
     299      meteor_ccn(:,:,:)=0.
    284300      rice(:,:) = 1.e-8
    285301      riceco2(:,:) = 1.e-11
    … …  
    320336            do ig=1,ngrid
    321337               masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
    322                ltemp1=abs(pression_meteo(:)-pplev(ig,l))
    323                ltemp2=abs(pression_meteo(:)-pplev(ig,l+1))
     338               ltemp1=abs(pression_meteor(:)-pplev(ig,l))
     339               ltemp2=abs(pression_meteor(:)-pplev(ig,l+1))
    324340               lebon1=minloc(ltemp1,DIM=1)
    325341               lebon2=minloc(ltemp2,DIM=1)
    … …  
    327343               mtemp(:)=0d0     !mtemp(100) : valeurs pour les 100bins
    328344               do ibin=1,100
    329                   mtemp(ibin)=sum(meteo(lebon1:lebon2,ibin))
     345                  mtemp(ibin)=sum(meteor(lebon1:lebon2,ibin))
    330346               enddo
    331                meteo_ccn(ig,l,:)=mtemp(:)/nelem/masse(ig,l) !Par kg air
     347               meteor_ccn(ig,l,:)=mtemp(:)/nelem/masse(ig,l) !Par kg air
    332348csi par m carre, x epaisseur/masse pour par kg/air
    333349               !write(*,*) "masse air ig l=",masse(ig,l)
    … …  
    388404                 m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
    389405                 n_aer(i) = n_aer(i) + 0.5 * No * n_derf +
    390      &                meteo_ccn(ig,l,i) !Ajout meteo_ccn particles
     406     &                meteor_ccn(ig,l,i) !Ajout meteor_ccn particles
    391407                 m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
    392408                 m_aer2(i)=4./3.*pi*rho_dust
    393      &                *meteo_ccn(ig,l,i)*rad_cldco2(i)*rad_cldco2(i)
     409     &                *meteor_ccn(ig,l,i)*rad_cldco2(i)*rad_cldco2(i)
    394410     &                *rad_cldco2(i)
    395411              enddo
    … …  
    599615     &       (zq(1:ngrid,1:nlay,igcm_dust_number) -
    600616     &       zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
    601         return
     617
    602618        end
    603619     

  • trunk/LMDZ.MARS/libf/phymars/massflowrateCO2.F

    r1816 r1884  
    55c     sublimation
    66c
    7 c     newton-raphson method
    8 c     CLASSICAL  (no SF etc.)
    9 c
    107c   inputs: Pressure (P), Temperature (T), saturation ratio (Sat),
    118c           particle radius (Radius), molecular mass of the atm (Matm)
    … …  
    1310c
    1411c   Authors: C. Listowski (2014) then J. Audouard (2016-2017)
     12c 
     13c
     14c   Updates:
     15c   --------
     16c   December 2017 - C. Listowski - Simplification of the derivation of
     17c   massflowrate by using explicit formula for surface temperature,
     18c   No Newton-Raphson routine anymore- see comment at relevant line
    1519c=======================================================================
    16       USE comcstfi_h
     20      USE comcstfi_h, ONLY: pi
    1721
    1822      implicit none
    1923
    2024      include "microphys.h"
    21 
    22 
    2325
    2426c   arguments: INPUT
    2527c   ----------
    26       REAL T,Matm
    27       REAL*8 SAT
    28       real P
    29       DOUBLE PRECISION Radius
     28      REAL,INTENT(in) :: T,Matm
     29      REAL*8,INTENT(in) :: SAT
     30      REAL,INTENT(in) :: P
     31      DOUBLE PRECISION,INTENT(in) :: Radius
    3032c   arguments: OUTPUT
    3133c   ----------
    32       DOUBLE PRECISION   Ic
     34      DOUBLE PRECISION,INTENT(out) ::   Ic
    3335c   Local Variables
    3436c   ----------
    35       DOUBLE PRECISION   Tcm
    36       DOUBLE PRECISION   T_inf, T_sup, T_dT
     37      DOUBLE PRECISION   Tsurf
    3738      DOUBLE PRECISION   C0,C1,C2
    3839      DOUBLE PRECISION   kmix,Lsub,cond
    39       DOUBLE PRECISION   rtsafe
    40       DOUBLE PRECISION   left, fval, dfval
    41 c    function for newton-raphson iterative method
    42 c    --------------------------
    43       EXTERNAL classical
    44          
    45       Tcm = 110!dble(T)    ! initialize pourquoi 0 et pas t(i)
    46       T_inf = 0d0
    47       T_sup = 800d0
    48       T_dT = 0.001  ! precision - mettre petit et limiter nb iteration?
    49      
    50 666   CONTINUE
    51 
    52       call  coefffunc(P,T,Sat,Radius,Matm,kmix,Lsub,C0,C1,C2)
    53 
    54       if (isnan(C0) .eqv. .true.)  C0=0d0   
    55 c     FIND SURFACE TEMPERATURE (Tc) : iteration sur t
     40      DOUBLE PRECISION   Ak
     41
     42      call  coefffunc(P,T,Sat,Radius,Matm,kmix,Lsub,C0,C1,C2,Ak)
     43   
     44      Tsurf = 1./C1*dlog(Sat/Ak) + T
     45
     46      !Note by CL - dec 2017 (see also technical note)
     47      !The above is a simplified version of Tsurf
     48      !compared to the one used by Listowski et al. 2014 (Ta), where a
     49      !Newton-Raphson routine must be used. Approximations made by
     50      !considering the orders of magnitude of the different factors lead to
     51      !simplification of the equation 5 of Listowski et al. (2014).
     52      !The error compared to the exact value determined by NR iterations
     53      !is less than 0.6% for all sizes, pressures, supersaturations
     54      !relevant to present Mars. Should also be ok for most conditions
     55      !in ancient Mars (However, needs to be double-cheked, as explained in
     56      !(Listowski et al. 2013,JGR)
     57
    5658      cond = 4.*pi*Radius*kmix
    57       Tcm = rtsafe(classical,T_inf,T_sup,T_dT,Radius,C0,C1,C2)
    58       if (Tcm.LE.0d0 .or. isnan(Tcm)) then ! unsignificant cases where S<<<Seq and Ncores <<1e-10
    59          Tcm = T
    60       endif
    61 c     THEN COMPUTE MASS FLUX Ic from FINAL Tsurface (Tc)
    62       Ic = (Tcm-T)
    63       Ic = cond*Ic/Lsub
    64 c regarder de combien varie la solution Ic entre Tcm et Tcm+T_dT
    65  
    66       RETURN
    67 
     59      Ic = cond*(Tsurf-T)/Lsub
     60 
    6861      END
    6962
    70 
    71 c****************************************************************
    72       FUNCTION rtsafe(funcd,x1,x2,xacc,Radius,C0,C1,C2)
    73 *
    74 *     Newton Raphsen routine (Numerical Recipe)
    75 *
    76 c****************************************************************
    77 
    78       implicit none
    79 
    80       INTEGER MAXIT
    81        DOUBLE PRECISION x1,x2,xacc
    82       DOUBLE PRECISION rtsafe
    83       DOUBLE PRECISION Radius
    84       DOUBLE PRECISION C0,C1,C2
    85 
    86 
    87       EXTERNAL funcd
    88 
    89       PARAMETER (MAXIT=10000)   !Maximum allowed number of iterations. Using a combination of Newton-Raphson and bisection,
    90                                 !find the root of a function bracketed between x1 and x2. The root, returned as the function value rtsafe,
    91                                 !will be refined until its accuracy is known within !!±xacc. funcd is a user-supplied subroutine which
    92                                 !returns both the function value and the first derivative of the function.
    93 
    94       INTEGER j
    95       DOUBLE PRECISION df,dx,dxold
    96       DOUBLE PRECISION f,fh,fl,temp,xh,xl
    97 
    98       call funcd(x1,fl,df,C0,C1,C2)
    99       call funcd(x2,fh,df,C0,C1,C2)
    100 
    101 
    102       if ((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.) ) then
    103          x1=0d0
    104          x2=1200d0         
    105          call funcd(x1,fl,df,C0,C1,C2)
    106          call funcd(x2,fh,df,C0,C1,C2)
    107 c         write(*,*) 'root must be bracketed in rtsafe'
    108       endif
    109 
    110 
    111       if (fl.eq.0.) then
    112          rtsafe=x1
    113          return
    114       else if (fh.eq.0.) then
    115          rtsafe=x2
    116          return
    117       else if (fl.lt.0.) then   !Orient the search so that f(xl) < 0.                                                                   
    118          xl=x1
    119          xh=x2
    120       else
    121         xh=x1
    122         xl=x2
    123       endif
    124       rtsafe = .5*(x1+x2)       !Initialize the guess for root,
    125       dxold  = abs(x2-x1)       !the stepsize before last,
    126       dx     = dxold            ! and the last step.
    127       call funcd(rtsafe,f,df,C0,C1,C2)
    128       DO 11 j=1,MAXIT           !Loop over allowed iterations.
    129 
    130          if (((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).gt.0.   ! Bisect if Newton out of range
    131      *    .or. abs(2.*f).gt.abs(dxold*df) ) then               ! or not decreasing fst enough
    132 
    133              dxold=dx
    134              dx=0.5*(xh-xl)
    135              rtsafe=xl+dx
    136              if (xl.eq.rtsafe) return                   !Change in root is negligible. Newton step acceptable. Take it.
    137          else
    138             dxold=dx
    139             dx=f/df
    140             temp=rtsafe
    141             rtsafe=rtsafe-dx
    142             if(temp.eq.rtsafe) return
    143          endif
    144 
    145         if(abs(dx).lt.xacc) return !Convergence criterion. The one new function evaluation per iteration. Maintain the bracket on the root.
    146          call funcd(rtsafe,f,df,C0,C1,C2)
    147         if(f.lt.0.) then
    148          xl=rtsafe
    149         else
    150          xh=rtsafe
    151         endif
    152 11    ENDDO
    153 
    154       rtsafe=0d0
    155       write(*,*) 'rtsafe exceeding maximum iterations,Tcm=',rtsafe
    156       return
    157 
    158       END
    159 
    160 
    16163c********************************************************************************
    162       subroutine classical(x,f,df,C0,C1,C2)
    163 c     Function to give as input to RTSAFE (NEWTON-RAPHOEN)
    164 c********************************************************************************
    165 
    166       implicit none
    167 
    168       DOUBLE PRECISION   x
    169       DOUBLE PRECISION  C0,C1,C2
    170       DOUBLE PRECISION f
    171       DOUBLE PRECISION  df
    172 
    173       f   = x + C0*exp(C1*x) - C2   ! start f
    174       df  = 1. + C0*C1*exp(C1*x)    ! start df
    175  
    176       return
    177       END 
    178 
    179 c********************************************************************************
    180 
    181       subroutine coefffunc(P,T,S,rc,Matm,kmix,Lsub,C0,C1,C2)
     64
     65      subroutine coefffunc(P,T,S,rc,Matm,kmix,Lsub,C0,C1,C2,Ak)
    18266
    18367c********************************************************************************
    18468c defini la fonction eq 6 papier 2014  (Listowski et al., 2014)
    18569      use tracer_mod, only: rho_ice_co2
    186       USE comcstfi_h
     70      USE comcstfi_h, ONLY: pi
    18771
    18872      implicit none
    … …  
    19175c   arguments: INPUT
    19276c   ----------------
    193       REAL P
    194       real T
    195       REAL*8 S
    196       double precision rc
    197       REAL Matm !g.mol-1 ( = mmean(ig,l) )
     77      REAL,INTENT(in) :: P
     78      REAL,INTENT(in) :: T
     79      REAL*8,INTENT(in) :: S
     80      DOUBLE PRECISION,INTENT(in) :: rc
     81      REAL,INTENT(in) :: Matm !g.mol-1 ( = mmean(ig,l) )
     82c   arguments: OUTPUT
     83c   ----------
     84      DOUBLE PRECISION,INTENT(out) ::  C0,C1,C2
     85      DOUBLE PRECISION,INTENT(out) ::  kmix,Lsub
    19886
    19987c   local:
    … …  
    20795      DOUBLE PRECISION vthatm,lpmt,rhoatm, vthco2 ! for Kn,th
    20896
    209 c   arguments: OUTPUT
    210 c   ----------
    211       DOUBLE PRECISION  C0,C1,C2
    212       DOUBLE PRECISION  kmix,Lsub
    21397
    21498c     DEFINE heat cap. J.kg-1.K-1 and To
    … …  
    265149      C1 = Lsub*mco2/(rgp*dble(T)**2)
    266150      C2 = dble(T) + dble(P)*mco2*Dv*Lsub*xinf/(kmix*rgp*dble(T))
    267       RETURN
     151
    268152      END
    269153
    … …  
    281165c     -----------
    282166     
    283       REAL P
    284       REAL Pbar                     !!! has to be in bar for the formula
    285       REAL T
     167      REAL,INTENT(in) :: P
     168      REAL,INTENT(in) :: T
    286169     
    287170c     output
    288171c     -----------
    289172
    290       DOUBLE PRECISION Diff
     173      DOUBLE PRECISION,INTENT(out) :: Diff
    291174   
    292175c      local
    293176c     -----------
    294177
     178      REAL Pbar                     !!! has to be in bar for the formula
    295179      DOUBLE PRECISION dva, dvb, Mab  ! Mab has to be in g.mol-1
    296180       
    … …  
    326210c     -----------
    327211         
    328          REAL T
    329          DOUBLE PRECISION x
    330          DOUBLE PRECISION rho !kg.m-3
     212         REAL,INTENT(in) :: T
     213         DOUBLE PRECISION,INTENT(in) :: x
     214         DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3
    331215
    332216c     outputs
    333217c     -----------
    334218
    335          DOUBLE PRECISION Kthmix
     219         DOUBLE PRECISION,INTENT(out) :: Kthmix
    336220
    337221c     local
    … …  
    393277        Kthmix = kco2*x1 /(x1*A11 + x2*A12) + kn2*x2 /(x1*A21 + x2*A22)
    394278        Kthmix = Kthmix*1e-3   ! from mW.m-1.K-1 to  W.m-1.K-1
    395         RETURN
    396279
    397280        END
    … …  
    412295c     -----------
    413296         
    414          REAL T
    415          DOUBLE PRECISION rho !kg.m-3
     297         REAL,INTENT(in) :: T
     298         DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3
    416299
    417300c     outputs
    418301c     -----------
    419302
    420          DOUBLE PRECISION kthn2
     303         DOUBLE PRECISION,INTENT(out) :: kthn2
    421304
    422305c     local
    … …  
    505388         kthn2 = k1 + k2
    506389
    507          RETURN
    508 
    509390         END
    510391
    … …  
    525406c     -----------
    526407         
    527       REAL T
     408      REAL,INTENT(in) :: T
    528409
    529410c     outputs
    530411c     -----------
    531412
    532       DOUBLE PRECISION visco
     413      DOUBLE PRECISION,INTENT(out) :: visco
    533414
    534415
    … …  
    585466c     -----------
    586467
    587       REAL T
    588       DOUBLE PRECISION rho
     468      REAL,INTENT(in) :: T
     469      DOUBLE PRECISION,INTENT(in) :: rho
    589470
    590471c     outputs
    591472c     -----------
    592473
    593       DOUBLE PRECISION kthco2
     474      DOUBLE PRECISION,INTENT(out) :: kthco2
    594475
    595476c     LOCAL
    … …  
    657538      kthco2 = (k1 + k2) *  Lambdac   ! mW
    658539
    659       RETURN
    660 
    661540      END
Note: See TracChangeset for help on using the changeset viewer.