Changeset 3464 for trunk/LMDZ.MARS


Ignore:
Timestamp:
Oct 18, 2024, 10:29:59 AM (5 weeks ago)
Author:
emillour
Message:

Mars PCM:
Some tidying in aeronomars:

  • make a jthermcalc_util.F to contain utility routines (used by jthermcal & jthermcalc_e107). Also add the possibility (turned off by default) in the interfast routine to do extra sanity checks.
  • turn chemthermos, euvheat, hrtherm, jthermcalc, jthermcalc_e107, paramphoto_compact and photochemistry into modules.

EM

Location:
trunk/LMDZ.MARS
Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/changelog.txt

    r3459 r3464  
    47094709- Addition in the "util" folder of two python scripts: one is to visualize easily any variable of a NetCDF file; the other is to display useful information about the variables of a NetCDF file to help for debugging.
    47104710- Move of the launching scripts for the 1D model from the "deftank" to the "util" folder.
     4711
     4712== 18/10/2024 == EM
     4713Some tidying in aeronomars:
     4714- make a jthermcalc_util.F to contain utility routines (used by jthermcal &
     4715  jthermcalc_e107). Also add the possibility (turned off by default) in the
     4716  interfast routine to do extra sanity checks.
     4717- turn chemthermos, euvheat, hrtherm, jthermcalc, jthermcalc_e107,
     4718  paramphoto_compact and photochemistry into modules

  • trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90

    r3461 r3464  
    3131      use conc_mod, only: mmean ! mean molecular mass of the atmosphere
    3232      use comcstfi_h, only: pi
     33      use chemthermos_mod, only: chemthermos
     34      use photochemistry_mod, only: photochemistry
    3335      use chemistrydata_mod, only: read_phototable
    3436      use photolysis_mod, only: init_photolysis, nphot

  • trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90

    r2615 r3464  
     1  MODULE chemthermos_mod
     2 
     3  IMPLICIT NONE
     4 
     5  CONTAINS
     6 
    17      SUBROUTINE chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp, &
    28           zdens,zpress,zlocal,zenit,ptimestep,zday,em_no,em_o2)
    … …  
    1218
    1319      use param_v4_h, only: Pno, Po2
    14       USE comcstfi_h
     20      use jthermcalc_e107_mod, only: jthermcalc_e107
     21      use paramfoto_compact_mod, only: paramfoto_compact
     22
    1523      IMPLICIT NONE
    1624!=======================================================================
    … …  
    3240!    ------------------
    3341!
    34 #include "callkeys.h"
     42      include "callkeys.h"
    3543!-----------------------------------------------------------------------
    3644!    Input/Output
    … …  
    540548      deallocate(rm)
    541549
    542       return
    543       end
    544 
    545 
    546 
    547 
     550      end subroutine chemthermos
     551
     552  END MODULE chemthermos_mod
     553
     554
     555

  • trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90

    r3158 r3464  
     1  MODULE euvheat_mod
     2 
     3  IMPLICIT NONE
     4 
     5  CONTAINS
     6
    17      SUBROUTINE euvheat(ngrid,nlayer,nq,pt,pdt,pplev,pplay,zzlay, &
    28           mu0,ptimestep,ptime,zday,pq,pdq,pdteuv)
    … …  
    713                            igcm_n, igcm_no, igcm_no2, igcm_n2d, mmol
    814      use conc_mod, only: rnew, cpnew
     15      use hrtherm_mod, only: hrtherm
    916      IMPLICIT NONE
    1017!=======================================================================
    … …  
    3138!    ------------------
    3239!
    33 #include "callkeys.h"
     40      include "callkeys.h"
    3441!-----------------------------------------------------------------------
    3542!    Input/Output
    … …  
    442449!      deallocate(rm)
    443450
    444       return
    445       end
     451      end subroutine euvheat
     452
     453  END MODULE euvheat_mod

  • trunk/LMDZ.MARS/libf/aeronomars/hrtherm.F

    r2615 r3464  
     1      MODULE hrtherm_mod
     2     
     3      IMPLICIT NONE
     4     
     5      CONTAINS
     6     
    17c**********************************************************************
    28
    … …  
    1117
    1218      use param_v4_h, only: ninter,nabs,jfotsout,fluxtop,freccen
     19      use jthermcalc_e107_mod, only: jthermcalc_e107
    1320
    1421      implicit none
    … …  
    132139      end do
    133140
    134       end
     141      end subroutine hrtherm
     142     
     143      END MODULE hrtherm_mod
    135144

  • trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F

    r2674 r3464  
     1      module jthermcalc_mod
     2     
     3      implicit none
     4     
     5      contains
     6
    17c**********************************************************************
    28
    … …  
    1622     .    co2crsc195,co2crsc295,t0,
    1723     .    jabsifotsintpar,ninter,nz2
    18 
     24      use jthermcalc_util, only: column, interfast
     25     
    1926      implicit none
    2027
    … …  
    10271034
    10281035
    1029       return
    1030 
    1031       end
    1032 
    1033 
    1034 
    1035 c**********************************************************************
    1036 c**********************************************************************
    1037 
    1038       subroutine column(ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit,
    1039      $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx,
    1040      $     n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
    1041 
    1042 c     nov 2002        fgg           first version
    1043 
    1044 c**********************************************************************
    1045 
    1046       use tracer_mod, only: igcm_o, igcm_co2, igcm_o2, igcm_h2,
    1047      &                      igcm_h2o_vap, igcm_h2o2, igcm_co, igcm_h,
    1048      &                      igcm_o3, igcm_n2, igcm_n, igcm_no, igcm_no2,
    1049      &                      mmol
    1050       use param_v4_h, only: radio,gg,masa,kboltzman,n_avog
    1051 
    1052       implicit none
    1053 
    1054 
    1055 c     common variables and constants
    1056       include 'callkeys.h'
    1057 
    1058 
    1059 
    1060 c    local parameters and variables
    1061 
    1062 
    1063 
    1064 c     input and output variables
    1065 
    1066       integer    ig,nlayer
    1067       integer    chemthermod
    1068       integer    nesptherm                      !# of species undergoing chemistry, input
    1069       real       rm(nlayer,nesptherm)         !densities (cm-3), input
    1070       real       tx(nlayer)                   !temperature profile, input
    1071       real       iz(nlayer+1)                 !height profile, input
    1072       real       zenit                          !SZA, input
    1073       real       co2colx(nlayer)              !column density of CO2 (cm^-2), output
    1074       real       o2colx(nlayer)               !column density of O2(cm^-2), output
    1075       real       o3pcolx(nlayer)              !column density of O(3P)(cm^-2), output
    1076       real       h2colx(nlayer)               !H2 column density (cm-2), output
    1077       real       h2ocolx(nlayer)              !H2O column density (cm-2), output
    1078       real       h2o2colx(nlayer)             !column density of H2O2(cm^-2), output
    1079       real       o3colx(nlayer)               !O3 column density (cm-2), output
    1080       real       n2colx(nlayer)               !N2 column density (cm-2), output
    1081       real       ncolx(nlayer)                !N column density (cm-2), output
    1082       real       nocolx(nlayer)               !NO column density (cm-2), output
    1083       real       cocolx(nlayer)               !CO column density (cm-2), output
    1084       real       hcolx(nlayer)                !H column density (cm-2), output
    1085       real       no2colx(nlayer)              !NO2 column density (cm-2), output
    1086 
    1087 
    1088 c     local variables
    1089 
    1090       real       xx
    1091       real       grav(nlayer)
    1092       real       Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2
    1093       real       Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2
    1094 
    1095       real       co2x(nlayer)
    1096       real       o2x(nlayer)
    1097       real       o3px(nlayer)
    1098       real       cox(nlayer)
    1099       real       hx(nlayer)
    1100       real       h2x(nlayer)
    1101       real       h2ox(nlayer)
    1102       real       h2o2x(nlayer)
    1103       real       o3x(nlayer)
    1104       real       n2x(nlayer)
    1105       real       nx(nlayer)
    1106       real       nox(nlayer)
    1107       real       no2x(nlayer)
    1108 
    1109       integer    i,j,k,icol,indexint          !indexes
    1110 
    1111 c     variables for optical path calculation
    1112 
    1113       integer    nz3
    1114 !      parameter  (nz3=nz*2)
    1115 
    1116       integer    jj
    1117       real*8      esp(nlayer*2)
    1118       real*8      ilayesp(nlayer*2)
    1119       real*8      szalayesp(nlayer*2)
    1120       integer     nlayesp
    1121       real*8      zmini
    1122       real*8      depth
    1123       real*8      espco2, espo2, espo3p, esph2, esph2o, esph2o2,espo3
    1124       real*8      espn2,espn,espno,espco,esph,espno2
    1125       real*8      rcmnz, rcmmini
    1126       real*8      szadeg
    1127 
    1128 
    1129       ! Tracer indexes in the thermospheric chemistry:
    1130       !!! ATTENTION. These values have to be identical to those in chemthermos.F90
    1131       !!! If the values are changed there, the same has to be done here  !!!
    1132       integer,parameter :: i_co2  =  1
    1133       integer,parameter :: i_co   =  2
    1134       integer,parameter :: i_o    =  3
    1135       integer,parameter :: i_o1d  =  4
    1136       integer,parameter :: i_o2   =  5
    1137       integer,parameter :: i_o3   =  6
    1138       integer,parameter :: i_h    =  7
    1139       integer,parameter :: i_h2   =  8
    1140       integer,parameter :: i_oh   =  9
    1141       integer,parameter :: i_ho2  = 10
    1142       integer,parameter :: i_h2o2 = 11
    1143       integer,parameter :: i_h2o  = 12
    1144       integer,parameter :: i_n    = 13
    1145       integer,parameter :: i_n2d  = 14
    1146       integer,parameter :: i_no   = 15
    1147       integer,parameter :: i_no2  = 16
    1148       integer,parameter :: i_n2   = 17
    1149 !      integer,parameter :: i_co2=1
    1150 !      integer,parameter :: i_o2=2
    1151 !      integer,parameter :: i_o=3
    1152 !      integer,parameter :: i_co=4
    1153 !      integer,parameter :: i_h=5
    1154 !      integer,parameter :: i_h2=8
    1155 !      integer,parameter :: i_h2o=9
    1156 !      integer,parameter :: i_h2o2=10
    1157 !      integer,parameter :: i_o3=12
    1158 !      integer,parameter :: i_n2=13
    1159 !      integer,parameter :: i_n=14
    1160 !      integer,parameter :: i_no=15
    1161 !      integer,parameter :: i_no2=17
    1162 
    1163 
    1164 c*************************PROGRAM STARTS*******************************
    1165 
    1166       nz3 = nlayer*2
    1167       do i=1,nlayer
    1168          xx = ( radio + iz(i) ) * 1.e5
    1169          grav(i) = gg * masa /(xx**2)
    1170       end do
    1171 
    1172       !Scale heights
    1173       xx = kboltzman * tx(nlayer) * n_avog / grav(nlayer)
    1174       Ho3p  = xx / mmol(igcm_o)
    1175       Hco2  = xx / mmol(igcm_co2)
    1176       Ho2   = xx / mmol(igcm_o2)
    1177       Hh2   = xx / mmol(igcm_h2)
    1178       Hh2o  = xx / mmol(igcm_h2o_vap)
    1179       Hh2o2 = xx / mmol(igcm_h2o2)
    1180       Hco   = xx / mmol(igcm_co)
    1181       Hh    = xx / mmol(igcm_h)
    1182       !Only if O3 chem. required
    1183       if(chemthermod.ge.1)
    1184      $     Ho3   = xx / mmol(igcm_o3)
    1185       !Only if N or ion chem.
    1186       if(chemthermod.ge.2) then
    1187          Hn2   = xx / mmol(igcm_n2)
    1188          Hn    = xx / mmol(igcm_n)
    1189          Hno   = xx / mmol(igcm_no)
    1190          Hno2  = xx / mmol(igcm_no2)
    1191       endif
    1192       ! first loop in altitude : initialisation
    1193       do i=nlayer,1,-1
    1194          !Column initialisation
    1195          co2colx(i)  = 0.
    1196          o2colx(i)   = 0.
    1197          o3pcolx(i)  = 0.
    1198          h2colx(i)   = 0.
    1199          h2ocolx(i)  = 0.
    1200          h2o2colx(i) = 0.
    1201          o3colx(i)   = 0.
    1202          n2colx(i)   = 0.
    1203          ncolx(i)    = 0.
    1204          nocolx(i)   = 0.
    1205          cocolx(i)   = 0.
    1206          hcolx(i)    = 0.
    1207          no2colx(i)  = 0.
    1208          !Densities
    1209          co2x(i)  = rm(i,i_co2)
    1210          o2x(i)   = rm(i,i_o2)
    1211          o3px(i)  = rm(i,i_o)
    1212          h2x(i)   = rm(i,i_h2)
    1213          h2ox(i)  = rm(i,i_h2o)
    1214          h2o2x(i) = rm(i,i_h2o2)
    1215          cox(i)   = rm(i,i_co)
    1216          hx(i)    = rm(i,i_h)
    1217          !Only if O3 chem. required
    1218          if(chemthermod.ge.1)
    1219      $        o3x(i)   = rm(i,i_o3)
    1220          !Only if Nitrogen of ion chemistry requested
    1221          if(chemthermod.ge.2) then
    1222             n2x(i)   = rm(i,i_n2)
    1223             nx(i)    = rm(i,i_n)
    1224             nox(i)   = rm(i,i_no)
    1225             no2x(i)  = rm(i,i_no2)
    1226          endif
    1227       enddo
    1228       ! second loop in altitude : column calculations
    1229       do i=nlayer,1,-1
    1230          !Routine to calculate the geometrical length of each layer
    1231          call espesor_optico_A(ig,i,nlayer,zenit,iz(i),nz3,iz,esp,
    1232      $         ilayesp,szalayesp,nlayesp, zmini)
    1233          if(ilayesp(nlayesp).eq.-1) then
    1234             co2colx(i)=1.e25
    1235             o2colx(i)=1.e25
    1236             o3pcolx(i)=1.e25
    1237             h2colx(i)=1.e25
    1238             h2ocolx(i)=1.e25
    1239             h2o2colx(i)=1.e25
    1240             o3colx(i)=1.e25
    1241             n2colx(i)=1.e25
    1242             ncolx(i)=1.e25
    1243             nocolx(i)=1.e25
    1244             cocolx(i)=1.e25
    1245             hcolx(i)=1.e25
    1246             no2colx(i)=1.e25
    1247          else
    1248             rcmnz = ( radio + iz(nlayer) ) * 1.e5
    1249             rcmmini = ( radio + zmini ) * 1.e5
    1250             !Column calculation taking into account the geometrical depth
    1251             !calculated before
    1252             do j=1,nlayesp
    1253                jj=ilayesp(j)
    1254                !Top layer
    1255                if(jj.eq.nlayer) then
    1256                   if(zenit.le.60.) then
    1257                      o3pcolx(i)=o3pcolx(i)+o3px(nlayer)*Ho3p*esp(j)
    1258      $                    *1.e-5
    1259                      co2colx(i)=co2colx(i)+co2x(nlayer)*Hco2*esp(j)
    1260      $                    *1.e-5
    1261                      h2o2colx(i)=h2o2colx(i)+
    1262      $                    h2o2x(nlayer)*Hh2o2*esp(j)*1.e-5
    1263                      o2colx(i)=o2colx(i)+o2x(nlayer)*Ho2*esp(j)
    1264      $                    *1.e-5
    1265                      h2colx(i)=h2colx(i)+h2x(nlayer)*Hh2*esp(j)
    1266      $                    *1.e-5
    1267                      h2ocolx(i)=h2ocolx(i)+h2ox(nlayer)*Hh2o*esp(j)
    1268      $                    *1.e-5                     
    1269                      cocolx(i)=cocolx(i)+cox(nlayer)*Hco*esp(j)
    1270      $                    *1.e-5
    1271                      hcolx(i)=hcolx(i)+hx(nlayer)*Hh*esp(j)
    1272      $                    *1.e-5
    1273                      !Only if O3 chemistry required
    1274                      if(chemthermod.ge.1) o3colx(i)=
    1275      $                    o3colx(i)+o3x(nlayer)*Ho3*esp(j)
    1276      $                    *1.e-5
    1277                      !Only if N or ion chemistry requested
    1278                      if(chemthermod.ge.2) then
    1279                         n2colx(i)=n2colx(i)+n2x(nlayer)*Hn2*esp(j)
    1280      $                    *1.e-5
    1281                         ncolx(i)=ncolx(i)+nx(nlayer)*Hn*esp(j)
    1282      $                       *1.e-5
    1283                         nocolx(i)=nocolx(i)+nox(nlayer)*Hno*esp(j)
    1284      $                       *1.e-5
    1285                         no2colx(i)=no2colx(i)+no2x(nlayer)*Hno2*esp(j)
    1286      $                       *1.e-5
    1287                      endif
    1288                   else if(zenit.gt.60.) then
    1289                      espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j)
    1290                      espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
    1291                      espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j)
    1292                      esph2  = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j)
    1293                      esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j)
    1294                      esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j)
    1295                      espco  = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j)
    1296                      esph   = sqrt((rcmnz+Hh)**2 -rcmmini**2)  - esp(j)
    1297                      !Only if O3 chemistry required
    1298                      if(chemthermod.ge.1)                     
    1299      $                   espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j)
    1300                      !Only if N or ion chemistry requested
    1301                      if(chemthermod.ge.2) then
    1302                         espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
    1303                         espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
    1304                         espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j)
    1305                         espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j)
    1306                      endif
    1307                      co2colx(i) = co2colx(i) + espco2*co2x(nlayer)
    1308                      o2colx(i)  = o2colx(i)  + espo2*o2x(nlayer)
    1309                      o3pcolx(i) = o3pcolx(i) + espo3p*o3px(nlayer)
    1310                      h2colx(i)  = h2colx(i)  + esph2*h2x(nlayer)
    1311                      h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(nlayer)
    1312                      h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(nlayer)
    1313                      cocolx(i)  = cocolx(i)  + espco*cox(nlayer)
    1314                      hcolx(i)   = hcolx(i)   + esph*hx(nlayer)
    1315                      !Only if O3 chemistry required
    1316                      if(chemthermod.ge.1)                     
    1317      $                  o3colx(i) = o3colx(i)  + espo3*o3x(nlayer)
    1318                      !Only if N or ion chemistry requested
    1319                      if(chemthermod.ge.2) then
    1320                         n2colx(i)  = n2colx(i)  + espn2*n2x(nlayer)
    1321                         ncolx(i)   = ncolx(i)   + espn*nx(nlayer)
    1322                         nocolx(i)  = nocolx(i)  + espno*nox(nlayer)
    1323                         no2colx(i) = no2colx(i) + espno2*no2x(nlayer)
    1324                      endif
    1325                   endif !Of if zenit.lt.60
    1326                !Other layers
    1327                else
    1328                   co2colx(i)  = co2colx(i) +
    1329      $                 esp(j) * (co2x(jj)+co2x(jj+1)) / 2.
    1330                   o2colx(i)   = o2colx(i) +
    1331      $                 esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
    1332                   o3pcolx(i)  = o3pcolx(i) +
    1333      $                 esp(j) * (o3px(jj)+o3px(jj+1)) / 2.
    1334                   h2colx(i)   = h2colx(i) +
    1335      $                 esp(j) * (h2x(jj)+h2x(jj+1)) / 2.
    1336                   h2ocolx(i)  = h2ocolx(i) +
    1337      $                 esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2.
    1338                   h2o2colx(i) = h2o2colx(i) +
    1339      $                 esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2.
    1340                   cocolx(i)   = cocolx(i) +
    1341      $                 esp(j) * (cox(jj)+cox(jj+1)) / 2.
    1342                   hcolx(i)    = hcolx(i) +
    1343      $                 esp(j) * (hx(jj)+hx(jj+1)) / 2.
    1344                   !Only if O3 chemistry required
    1345                   if(chemthermod.ge.1)
    1346      $                 o3colx(i) = o3colx(i) +
    1347      $                 esp(j) * (o3x(jj)+o3x(jj+1)) / 2.
    1348                   !Only if N or ion chemistry requested
    1349                   if(chemthermod.ge.2) then
    1350                      n2colx(i)   = n2colx(i) +
    1351      $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
    1352                      ncolx(i)    = ncolx(i) +
    1353      $                    esp(j) * (nx(jj)+nx(jj+1)) / 2.
    1354                      nocolx(i)   = nocolx(i) +
    1355      $                    esp(j) * (nox(jj)+nox(jj+1)) / 2.
    1356                      no2colx(i)  = no2colx(i) +
    1357      $                    esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
    1358                   endif
    1359                endif  !Of if jj.eq.nlayer
    1360             end do    !Of do j=1,nlayesp
    1361          endif        !Of ilayesp(nlayesp).eq.-1
    1362       enddo           !Of do i=nlayer,1,-1
    1363 
    1364       return
    1365 
    1366 
    1367       end
    1368 
    1369 
    1370 c**********************************************************************
    1371 c**********************************************************************
    1372 
    1373       subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup)
    1374 C
    1375 C subroutine to perform linear interpolation in pressure from 1D profile
    1376 C escin(nl) sampled on pressure grid pin(nl) to profile
    1377 C escout(nlayer) on pressure grid p(nlayer).
    1378 C
    1379       real*8,intent(out) :: wm(nlayer),wp(nlayer) ! interpolation weights
    1380       integer,intent(out) :: nm(nlayer) ! index of nearest point
    1381       real*8,intent(in) :: pin(nl),p(nlayer)
    1382       real*8,intent(in) :: limup,limdown
    1383       integer,intent(in) :: nl,nlayer
    1384       integer :: n1,n,np,nini
    1385       nini=1
    1386       do n1=1,nlayer
    1387          if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then
    1388             wm(n1) = 0.d0
    1389             wp(n1) = 0.d0
    1390          else
    1391             do n = nini,nl-1
    1392                if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
    1393                   nm(n1)=n
    1394                   np=n+1
    1395                   wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
    1396                   wp(n1)=1.d0 - wm(n1)
    1397                   nini = n
    1398                   exit
    1399                endif
    1400             enddo
    1401          endif
    1402       enddo
    1403 
    1404       end
    1405 
    1406 
    1407 c**********************************************************************
    1408 c**********************************************************************
    1409 
    1410       subroutine espesor_optico_A (ig,capa,nlayer, szadeg,z,
    1411      @                   nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
    1412 
    1413 c     fgg              nov 03      Adaptation to Martian model
    1414 c     malv             jul 03      Corrected z grid. Split in alt & frec codes
    1415 c     fgg              feb 03      first version
    1416 *************************************************************************
    1417 
    1418       use param_v4_h, only: radio
    1419       implicit none
    1420 
    1421 c     arguments
    1422 
    1423       real        szadeg                ! I. SZA [rad]
    1424       real        z                     ! I. altitude of interest [km]
    1425       integer     nz3,ig,nlayer         ! I. dimension of esp, ylayesp, etc...
    1426                                         !  (=2*nlayer= max# of layers in ray path)
    1427       real     iz(nlayer+1)              ! I. Altitude of each layer
    1428       real*8        esp(nz3)            ! O. layer widths after geometrically
    1429                                         !    amplified; in [cm] except at TOA
    1430                                         !    where an auxiliary value is used
    1431       real*8        ilayesp(nz3)        ! O. Indexes of layers along ray path
    1432       real*8        szalayesp(nz3)      ! O. SZA [deg]    "     "       "
    1433       integer       nlayesp
    1434 !      real*8        nlayesp             ! O. # layers along ray path at this z
    1435       real*8        zmini               ! O. Minimum altitud of ray path [km]
    1436 
    1437 
    1438 c     local variables and constants
    1439 
    1440         integer     j,i,capa
    1441         integer     jmin                  ! index of min.altitude along ray path
    1442         real*8      szarad                ! SZA [deg]
    1443         real*8      zz
    1444         real*8      diz(nlayer+1)
    1445         real*8      rkmnz                 ! distance TOA to center of Planet [km]
    1446         real*8      rkmmini               ! distance zmini to center of P [km]
    1447         real*8      rkmj                  ! intermediate distance to C of P [km]
    1448 c external function
    1449         external  grid_R8          ! Returns index of layer containing the altitude
    1450                                 ! of interest, z; for example, if
    1451                                 ! zkm(i)=z or zkm(i)<z<zkm(i+1) => grid(z)=i
    1452         integer   grid_R8
    1453 
    1454 *************************************************************************     
    1455         szarad = dble(szadeg)*3.141592d0/180.d0
    1456         zz=dble(z)
    1457         do i=1,nlayer
    1458            diz(i)=dble(iz(i))
    1459         enddo
    1460         do j=1,nz3
    1461           esp(j) = 0.d0
    1462           szalayesp(j) = 777.d0
    1463           ilayesp(j) = 0
    1464         enddo
    1465         nlayesp = 0
    1466 
    1467         ! First case: szadeg<60
    1468         ! The optical thickness will be given by  1/cos(sza)
    1469         ! We deal with 2 different regions:
    1470         !   1: First, all layers between z and ztop ("upper part of ray")
    1471         !   2: Second, the layer at ztop
    1472         if(szadeg.lt.60.d0) then
    1473 
    1474            zmini = zz
    1475            if(abs(zz-diz(nlayer)).lt.1.d-3) goto 1357
    1476            ! 1st Zone: Upper part of ray
    1477            !
    1478            do j=grid_R8(zz,diz,nlayer),nlayer-1
    1479              nlayesp = nlayesp + 1
    1480              ilayesp(nlayesp) = j
    1481              esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad)        ! [km]
    1482              esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
    1483              szalayesp(nlayesp) = szadeg
    1484            end do
    1485 
    1486            !
    1487            ! 2nd Zone: Top layer
    1488  1357      continue
    1489            nlayesp = nlayesp + 1
    1490            ilayesp(nlayesp) = nlayer
    1491            esp(nlayesp) = 1.d0 / cos(szarad)         ! aux. non-dimens. factor
    1492            szalayesp(nlayesp) = szadeg
    1493 
    1494 
    1495         ! Second case:  60 < szadeg < 90
    1496         ! The optical thickness is evaluated.
    1497         !    (the magnitude of the effect of not using cos(sza) is 3.e-5
    1498         !     for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately)
    1499         ! We deal with 2 different regions:
    1500         !   1: First, all layers between z and ztop ("upper part of ray")
    1501         !   2: Second, the layer at ztop ("uppermost layer")
    1502         else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
    1503 
    1504            zmini=(radio+zz)*sin(szarad)-radio
    1505            rkmmini = radio + zmini
    1506 
    1507            if(abs(zz-diz(nlayer)).lt.1.d-4) goto 1470
    1508 
    1509            ! 1st Zone: Upper part of ray
    1510            !
    1511            do j=grid_R8(zz,diz,nlayer),nlayer-1
    1512               nlayesp = nlayesp + 1
    1513               ilayesp(nlayesp) = j
    1514               esp(nlayesp) =
    1515      &             sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
    1516      &             sqrt( (radio+diz(j))**2 - rkmmini**2 )           ! [km]
    1517               esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
    1518               rkmj = radio+diz(j)
    1519               szalayesp(nlayesp) = asin( rkmmini/rkmj )             ! [rad]
    1520               szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg]
    1521            end do
    1522  1470      continue
    1523            ! 2nd Zone:  Uppermost layer of ray.
    1524            !
    1525            nlayesp = nlayesp + 1
    1526            ilayesp(nlayesp) = nlayer
    1527            rkmnz = radio+diz(nlayer)
    1528            esp(nlayesp)  =  sqrt( rkmnz**2 - rkmmini**2 )       ! aux.factor[km]
    1529            esp(nlayesp)  =  esp(nlayesp) * 1.d5                 ! aux.f. [cm]
    1530            szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
    1531            szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
    1532 
    1533 
    1534         ! Third case:  szadeg > 90
    1535         ! The optical thickness is evaluated.
    1536         ! We deal with 5 different regions:
    1537         !   1: all layers between z and ztop ("upper part of ray")
    1538         !   2: the layer at ztop ("uppermost layer")
    1539         !   3: the lowest layer, at zmini
    1540         !   4: the layers increasing from zmini to z (here SZA<90)
    1541         !   5: the layers decreasing from z to zmini (here SZA>90)
    1542         else if(szadeg.gt.90.d0) then
    1543 
    1544            zmini=(radio+zz)*sin(szarad)-radio
    1545            !zmini should be lower than zz, as SZA<90. However, in situations
    1546            !where SZA is very close to 90, rounding errors can make zmini
    1547            !slightly higher than zz, causing problems in the determination
    1548            !of the jmin index. A correction is implemented in the determination
    1549            !of jmin, some lines below
    1550            rkmmini = radio + zmini
    1551 
    1552            if(zmini.lt.diz(1)) then         ! Can see the sun?  No => esp(j)=inft
    1553              nlayesp = nlayesp + 1
    1554              ilayesp(nlayesp) = - 1     ! Value to mark "no sun on view"
    1555 !             esp(nlayesp) = 1.e30
    1556 
    1557            else
    1558               jmin=grid_R8(zmini,diz,nlayer)+1
    1559               !Correction for possible rounding errors when SZA very close
    1560               !to 90 degrees
    1561               if(jmin.gt.grid_R8(zz,diz,nlayer)) then
    1562                  write(*,*)'jthermcalc warning: possible rounding error'
    1563                  write(*,*)'point,sza,layer:',ig,szadeg,capa
    1564                  jmin=grid_R8(zz,diz,nlayer)
    1565               endif
    1566 
    1567               if(abs(zz-diz(nlayer)).lt.1.d-4) goto 9876
    1568 
    1569               ! 1st Zone: Upper part of ray
    1570               !
    1571               do j=grid_R8(zz,diz,nlayer),nlayer-1
    1572                 nlayesp = nlayesp + 1
    1573                 ilayesp(nlayesp) = j
    1574                 esp(nlayesp) =
    1575      $                sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
    1576      $                sqrt( (radio+diz(j))**2 - rkmmini**2 )          ! [km]
    1577                 esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
    1578                 rkmj = radio+diz(j)
    1579                 szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
    1580                 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592      ! [deg]
    1581               end do
    1582 
    1583  9876         continue
    1584               ! 2nd Zone:  Uppermost layer of ray.
    1585               !
    1586               nlayesp = nlayesp + 1
    1587               ilayesp(nlayesp) = nlayer
    1588               rkmnz = radio+diz(nlayer)
    1589               esp(nlayesp) =  sqrt( rkmnz**2 - rkmmini**2 )      !aux.factor[km]
    1590               esp(nlayesp) = esp(nlayesp) * 1.d5                 !aux.f.[cm]
    1591               szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
    1592               szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
    1593 
    1594               ! 3er Zone: Lowestmost layer of ray
    1595               !
    1596               if ( jmin .ge. 2 ) then      ! If above the planet's surface
    1597                 j=jmin-1
    1598                 nlayesp = nlayesp + 1
    1599                 ilayesp(nlayesp) = j
    1600                 esp(nlayesp) = 2. *
    1601      $                 sqrt( (radio+diz(j+1))**2 -rkmmini**2 )       ! [km]
    1602                 esp(nlayesp) = esp(nlayesp) * 1.d5                   ! [cm]
    1603                 rkmj = radio+diz(j+1)
    1604                 szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
    1605                 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
    1606               endif
    1607 
    1608               ! 4th zone: Lower part of ray, increasing from zmin to z
    1609               !    ( layers with SZA < 90 deg )
    1610               do j=jmin,grid_R8(zz,diz,nlayer)-1
    1611                 nlayesp = nlayesp + 1
    1612                 ilayesp(nlayesp) = j
    1613                 esp(nlayesp) =
    1614      $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
    1615      $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )       ! [km]
    1616                 esp(nlayesp) = esp(nlayesp) * 1.d5                     ! [cm]
    1617                 rkmj = radio+diz(j)
    1618                 szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
    1619                 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
    1620               end do
    1621 
    1622               ! 5th zone: Lower part of ray, decreasing from z to zmin
    1623               !    ( layers with SZA > 90 deg )
    1624               do j=grid_R8(zz,diz,nlayer)-1, jmin, -1
    1625                 nlayesp = nlayesp + 1
    1626                 ilayesp(nlayesp) = j
    1627                 esp(nlayesp) =
    1628      $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
    1629      $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )        ! [km]
    1630                 esp(nlayesp) = esp(nlayesp) * 1.d5                      ! [cm]
    1631                 rkmj = radio+diz(j)
    1632                 szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )          ! [rad]
    1633                 szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg]
    1634               end do
    1635 
    1636            end if
    1637 
    1638         end if
    1639 
    1640 
    1641         return
    1642 
    1643         end
    1644 
    1645 
    1646 
    1647 c**********************************************************************
    1648 c***********************************************************************
    1649 
    1650         function grid_R8 (z, zgrid, nz)
    1651 
    1652 c Returns the index where z is located within vector zgrid
    1653 c The vector zgrid must be monotonously increasing, otherwise program stops.
    1654 c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops.
    1655 c
    1656 c FGG     Aug-2004     Correct z.lt.zgrid(i) to .le.
    1657 c MALV    Jul-2003
    1658 c***********************************************************************
    1659 
    1660         implicit none
    1661 
    1662 c Arguments
    1663         integer   nz
    1664         real*8      z
    1665         real*8      zgrid(nz)
    1666         integer   grid_R8
    1667 
    1668 c Local 
    1669         integer   i, nz1, nznew
    1670 
    1671 c*** CODE START
    1672 
    1673         if ( z .lt. zgrid(1)  ) then
    1674            write (*,*) ' GRID/ z outside bounds of zgrid '
    1675            write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
    1676            z = zgrid(1)
    1677            write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)'
    1678            write(*,*) 'Please check values of z and zgrid above'
    1679         endif
    1680         if (z .gt. zgrid(nz) ) then
    1681            write (*,*) ' GRID/ z outside bounds of zgrid '
    1682            write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
    1683            z = zgrid(nz)
    1684            write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)'
    1685            write(*,*) 'Please check values of z and zgrid above'
    1686         endif
    1687         if ( nz .lt. 2 ) then
    1688            write (*,*) ' GRID/ zgrid needs 2 points at least ! '
    1689            stop ' Serious error in GRID.F '
    1690         endif
    1691         if ( zgrid(1) .ge. zgrid(nz) ) then
    1692            write (*,*) ' GRID/ zgrid must increase with index'
    1693            stop ' Serious error in GRID.F '
    1694         endif
    1695 
    1696         nz1 = 1
    1697         nznew = nz/2
    1698         if ( z .gt. zgrid(nznew) ) then
    1699            nz1 = nznew
    1700            nznew = nz
    1701         endif
    1702         do i=nz1+1,nznew
    1703            if ( z. eq. zgrid(i) ) then
    1704               grid_R8=i
    1705               return
    1706               elseif ( z .le. zgrid(i) ) then
    1707               grid_R8 = i-1
    1708               return
    1709            endif
    1710         enddo
    1711         grid_R8 = nz
    1712         return
    1713 
    1714 
    1715 
    1716         end
    1717 
    1718 
    1719 
    1720 !c***************************************************
    1721 !c***************************************************
    1722 
    1723       subroutine flujo(date)
    1724 
    1725 
    1726 !c     fgg           nov 2002     first version
    1727 !c***************************************************
    1728 
    1729       use comsaison_h, only: dist_sol
    1730       use param_v4_h, only: ninter,
    1731      .                      fluxtop, ct1, ct2, p1, p2
    1732       implicit none
    1733 
    1734 
    1735 !     common variables and constants
    1736       include "callkeys.h"
    1737 
    1738 
    1739 !     Arguments
    1740 
    1741       real date
    1742 
    1743 
    1744 !     Local variable and constants
    1745 
    1746       integer i
    1747       integer inter
    1748       real    nada
    1749 
    1750 !c*************************************************
    1751 
    1752       if(date.lt.1985.) date=1985.
    1753       if(date.gt.2001.) date=2001.
     1036      end subroutine jthermcalc
    17541037     
    1755       do i=1,ninter
    1756          fluxtop(i)=1.
    1757          !Variation of solar flux with 11 years solar cycle
    1758          !For more details, see Gonzalez-Galindo et al. 2005
    1759          !To be improved in next versions
    1760         if(i.le.24.and.solvarmod.eq.0) then
    1761             fluxtop(i)=(((ct1(i)+p1(i)*date)/2.)                 
    1762      $           *sin(2.*3.1416/11.*(date-1985.-3.1416))         
    1763      $           +(ct2(i)+p2(i)*date)+1.)*fluxtop(i)
    1764          end if
    1765          fluxtop(i)=fluxtop(i)*(1.52/dist_sol)**2
    1766       end do
    1767      
    1768       return
    1769       end
     1038      end module jthermcalc_mod
     1039
     1040
     1041

  • trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F

    r2808 r3464  
     1      module jthermcalc_e107_mod
     2     
     3      implicit none
     4     
     5      contains
     6
    17c**********************************************************************
    28
    … …  
    1925     .    coefit0,coefit1,coefit2,coefit3,coefit4
    2026      use comsaison_h, only: dist_sol
     27      use jthermcalc_util, only: column, interfast
    2128
    2229      implicit none
    … …  
    11321139      jfotsout(:,:,:)=jfotsout(:,:,:)*(1.52/dist_sol)**2
    11331140
    1134       end
    1135 
    1136 
    1137 
    1138 
    1139 
    1140 
     1141      end subroutine jthermcalc_e107
     1142     
     1143      end module jthermcalc_e107_mod
     1144
     1145
     1146
     1147
     1148
     1149

  • trunk/LMDZ.MARS/libf/aeronomars/paramfoto_compact.F

    r2615 r3464  
     1      MODULE paramfoto_compact_mod
     2     
     3      IMPLICIT NONE
     4     
     5      CONTAINS
     6
    17c**********************************************************************
    28
    … …  
    639645cccccc End altitude loop
    640646
    641       return
    642 
    643 
    644       end
     647      end subroutine paramfoto_compact
    645648
    646649
    … …  
    693696         if ( c_output.lt.1.d-30) c_output=1.d-30
    694697
    695 
    696          return
    697698c END
    698          end
     699         end subroutine implicito
    699700
    700701
    … …  
    772773      return                                         
    773774
    774       end
     775      end function ionsec_nplus
    775776
    776777
    … …  
    840841      return                                         
    841842
    842       end  
     843      end function ionsec_n2plus
    843844
    844845
    … …  
    915916      return                                         
    916917
    917       end
     918      end function ionsec_oplus
    918919
    919920
    … …  
    993994      return                                         
    994995
    995       end  
     996      end function ionsec_coplus
    996997     
    997998
    … …  
    10691070      return                                         
    10701071
    1071       end  
     1072      end function ionsec_co2plus
    10721073     
    10731074     
    … …  
    11391140      return                                         
    11401141
    1141       end  
     1142      end function ionsec_o2plus
    11421143
    11431144
    … …  
    13461347
    13471348
    1348       return
    1349      
    1350 
    1351       end
     1349      end subroutine phdisrate
    13521350
    13531351
    … …  
    19891987        endif    !Of if(chemthermod.eq.3)
    19901988
    1991         return
    1992         end
     1989        end subroutine getch
    19931990
    19941991
    … …  
    20502047
    20512048
    2052       external  ionsec_nplus
    2053       real*8    ionsec_nplus
    2054 
    2055       external  ionsec_n2plus
    2056       real*8    ionsec_n2plus
    2057 
    2058       external  ionsec_oplus
    2059       real*8    ionsec_oplus
     2049!      external  ionsec_nplus
     2050!      real*8    ionsec_nplus
     2051
     2052!      external  ionsec_n2plus
     2053!      real*8    ionsec_n2plus
     2054
     2055!      external  ionsec_oplus
     2056!      real*8    ionsec_oplus
    20602057     
    2061       external  ionsec_coplus
    2062       real*8    ionsec_coplus
    2063 
    2064       external  ionsec_co2plus
    2065       real*8    ionsec_co2plus
    2066 
    2067       external  ionsec_o2plus
    2068       real*8    ionsec_o2plus
     2058!      external  ionsec_coplus
     2059!      real*8    ionsec_coplus
     2060
     2061!      external  ionsec_co2plus
     2062!      real*8    ionsec_co2plus
     2063
     2064!      external  ionsec_o2plus
     2065!      real*8    ionsec_o2plus
    20692066
    20702067
    … …  
    26882685         endif  !Of chemthermod.eq.3
    26892686
    2690          return
    26912687c END
    2692          end
     2688         end subroutine lifetimes
    26932689
    26942690
    … …  
    28682864      deltat = tminaux
    28692865
    2870       return
    28712866c END
    2872       end
     2867      end subroutine timemarching
    28732868
    28742869
    … …  
    29252920      integer   j
    29262921
    2927       external  ionsec_nplus
    2928       real*8    ionsec_nplus
    2929 
    2930       external  ionsec_n2plus
    2931       real*8    ionsec_n2plus
    2932 
    2933       external  ionsec_oplus
    2934       real*8    ionsec_oplus
    2935      
    2936       external  ionsec_coplus
    2937       real*8    ionsec_coplus
    2938 
    2939       external  ionsec_co2plus
    2940       real*8    ionsec_co2plus
    2941 
    2942       external  ionsec_o2plus
    2943       real*8    ionsec_o2plus
    2944 
    2945       logical firstcall
    2946       save firstcall
    2947 
     2922      logical,save :: firstcall=.true.
    29482923!$OMP THREADPRIVATE(firstcall)
    2949 
    2950       data firstcall /.true./
    29512924
    29522925ccccccccccccccc CODE STARTS
    … …  
    41004073      endif    !Of chemthermod.eq.3
    41014074
    4102       return
    41034075c END
    4104       end
     4076      end subroutine prodsandlosses
    41054077
    41064078
    … …  
    42434215      parameter (cociopt=1)
    42444216
    4245 c external functions
    4246 c
    4247       external cociente
    4248       real*8   cociente
    4249 
    4250       external ionsec_nplus
    4251       real*8   ionsec_nplus
    4252 
    4253       external ionsec_n2plus
    4254       real*8   ionsec_n2plus
    4255 
    4256       external ionsec_oplus
    4257       real*8   ionsec_oplus
    4258 
    4259       external ionsec_coplus
    4260       real*8   ionsec_coplus
    4261 
    4262       external ionsec_co2plus
    4263       real*8   ionsec_co2plus
    4264 
    4265       external ionsec_o2plus
    4266       real*8   ionsec_o2plus
    4267 
    4268       external avg
    4269       real*8   avg
    4270 
    4271       external dif
    4272       real*8   dif
    42734217
    42744218      real*8 log1
    … …  
    50835027
    50845028
    5085       return
    5086       end
     5029      end subroutine EF_oscilacion
    50875030
    50885031!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    50935036        avg = (log1+log2+log3)*0.333
    50945037        return
    5095         end
     5038        end function avg
    50965039!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    50975040        function dif(log1,log2,log3,avg)
    … …  
    51045047     &                abs(log3-avg) ) * 0.333
    51055048        return
    5106         end
     5049        end function dif
    51075050!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    51085051
    … …  
    51775120        return                                         
    51785121c END
    5179         end
     5122        end function cociente
     5123
     5124      END MODULE paramfoto_compact_mod

  • trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90

    r3462 r3464  
     1module photochemistry_mod
     2
     3implicit none
     4
     5contains
     6
    17!****************************************************************
    28!
    … …  
    2329
    2430use param_v4_h, only: jion
     31use jthermcalc_e107_mod, only: jthermcalc_e107
     32use paramfoto_compact_mod, only: phdisrate
    2533
    2634implicit none
    … …  
    43394347
    43404348end subroutine photochemistry
     4349
     4350end module photochemistry_mod

  • trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F

    r3443 r3464  
    1818      use moldiff_red_mod, only: moldiff_red ! old molecular diffusion scheme
    1919      use moldiff_MPF_mod, only: moldiff_MPF ! new molecular diffusion scheme
     20      use euvheat_mod, only: euvheat
    2021      use conc_mod, only: rnew, cpnew
    2122      USE comcstfi_h, only: r, cpp
Note: See TracChangeset for help on using the changeset viewer.