Ignore:
Timestamp:
Dec 4, 2024, 4:04:54 PM (3 weeks ago)
Author:
jbclement
Message:

PEM:
Removing unecessary module/subroutine "interpol_TI_PEM2PCM.F90" + Few small corrections/cleanings.
JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
1 deleted
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90

    r3527 r3532  
    1111!!!
    1212!!!
    13 !!! Author: EV, updated NS MSIM dynamical program for the PEM 
     13!!! Author: EV, updated NS MSIM dynamical program for the PEM
    1414!!!
    1515!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    3737  real(8)  ssi_depth_in, ssi_depth, T1
    3838  real(8), dimension(NP) :: zdepthF, zdepthE, zdepthT, zdepthG
    39   real(8), dimension(NMAX,NP) :: porefill, porefill_in 
     39  real(8), dimension(NMAX,NP) :: porefill, porefill_in
    4040  real(8), dimension(nz) :: Tb
    4141  real(8), dimension(NP) :: Tmean1, Tmean3, avrho1
    … …  
    6868     !read(21,*) latitude(k),albedo(k),thIn(k),htopo(k)
    6969     ! empirical relation from Mellon & Jakosky
    70      rhoc(k) = 800.*(150.+100.*sqrt(34.2+0.714*thIn(k))) 
     70     rhoc(k) = 800.*(150.+100.*sqrt(34.2+0.714*thIn(k)))
    7171  enddo
    7272  !close(21)
    … …  
    8080    z(iloop + 1) = l1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
    8181  enddo
    82  
     82
    8383
    8484  !open(unit=30,file='z.'//ext,action='write',status='unknown')
    … …  
    9595  timestep = 1 ! must be integer fraction of 1 ka
    9696  icetime = -tmax-timestep  ! earth years
    97  
    98   ! initializations 
     97
     98  ! initializations
    9999  !Tb = -9999.
    100100  zdepthF(:) = -9999.
    … …  
    152152  !print *,'Zt0=  ',ZdepthT
    153153     call icelayer_mars(timestep,nz,NP,thIn,rhoc,z,porosity,pfrost,Tb,zdepthF, &
    154           & zdepthE,porefill(1:nz,:),Tmean1,Tmean3,zdepthG, & 
     154          & zdepthE,porefill(1:nz,:),Tmean1,Tmean3,zdepthG, &
    155155          & albedo,p0,icefrac,zdepthT,avrho1, &
    156156          & avrho1prescribed)
    … …  
    164164         !write(36,*) icetime,latitude(k),zdepthF(k),porefill(1:nz,k)
    165165           ! compact output format
    166  !          write(36,'(f10.0,2x,f7.3,1x,f11.5,1x)',advance='no') & 
     166 !          write(36,'(f10.0,2x,f7.3,1x,f11.5,1x)',advance='no') &
    167167 !               & icetime,latitude(k),zdepthF(k)
    168168 !          call compactoutput(36,porefill(:,k),nz)
    … …  
    179179
    180180!501 format (f10.0,2x,f7.3,2x,f10.4,2(2x,f6.2),2x,f9.3,2x,g11.4)
    181  
     181
    182182end subroutine dyn_ss_ice_m
    183183
    … …  
    232232     endif
    233233     newti = sqrt(newrhoc*k)
    234      
     234
    235235  case (2)  ! massive ice (pure or dirty ice)
    236236     newrhoc = rhocobs*(1.-icefrac)/(1.-porosity) + icefrac*icedensity*cice
    237237     k = icefrac*kice + (1.-icefrac)*kw
    238238     newti = sqrt(newrhoc*k)
    239  
     239
    240240  case (3)  ! all ice, special case of layertype 2, which doesn't use porosity
    241241     newrhoc = icedensity*cice
    242      k = kice 
    243      newti = sqrt(newrhoc*k)
    244  
     242     k = kice
     243     newti = sqrt(newrhoc*k)
     244
    245245  case (4)  ! pores completely filled with ice, special case of layertype 1
    246246     newrhoc = rhocobs + porosity*icedensity*cice
    247      k = porosity*kice + kobs ! option A, end-member case of type 1, option A 
     247     k = porosity*kice + kobs ! option A, end-member case of type 1, option A
    248248     !k = kw*kice/((1-porosity)*kice+porosity*kw) ! option B, harmonic average
    249249     newti = sqrt(newrhoc*k)
    … …  
    257257  case default
    258258     error stop 'invalid layer type'
    259      
     259
    260260  end select
    261  
     261
    262262end subroutine soilthprop
    263263
    264264
    265265!=======================================================================
    266      
     266
    267267      real*8 function frostpoint(p)
    268268!     inverse of psv
    … …  
    271271      implicit none
    272272      real*8 p
    273      
     273
    274274!-----inverse of parametrization 1
    275275!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt
    … …  
    278278!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
    279279!      frostpoint = 1./(1./Tt-R/DHsub*log(p/pt))
    280      
     280
    281281!-----inverse of parametrization 2
    282282!     inverse of eq. (2) in Murphy & Koop (2005)
    … …  
    288288!     eq. (8) in Murphy & Koop (2005)
    289289!      frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p))
    290      
     290
    291291      end function frostpoint
    292292

  • trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90

    r3527 r3532  
    1111  real(8), parameter :: emiss0 = 1.  ! emissivity of dry surface
    1212  integer, parameter :: EQUILTIME =15 ! [Mars years]
    13  
     13
    1414  integer, parameter :: NMAX = 1000
    1515
    … …  
    5959     B = Diff*bigstep*86400.*365.24/(porosity*927.)
    6060     !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(Tb(),'h2o'))
    61      
     61
    6262     typeT = -9
    6363     if (zdepthT(k)>=0. .and. zdepthT(k)<z(nz)) then
    … …  
    7272   !  call assignthermalproperties(nz,thIn(k),rhoc(k),ti,rhocv, &
    7373   !       & typeT,icefrac,porosity,porefill(:,k))
    74      
     74
    7575     !----run thermal model
    76      call ajsub_mars(typeT, albedo(k), pfrost(k), nz, z, & 
     76     call ajsub_mars(typeT, albedo(k), pfrost(k), nz, z, &
    7777          &     ti, rhocv, fracIR, fracDust, p0(k), &
    7878          &     avdrho(k), avdrhoP(k), avrho1(k), Tb(k), zdepthE(k), typeF, &
    … …  
    9090        jump = -9
    9191     endif
    92      if (zdepthT(k)>=0. .and. avdrho(k)*avdrho_old(k)<0.) then 
     92     if (zdepthT(k)>=0. .and. avdrho(k)*avdrho_old(k)<0.) then
    9393        write(34,*) '# zdepth arrested'
    9494        if (jump>1) write(34,*) '# previous step was too large',jump
    … …  
    9999     avdrho_old(k) = avdrho(k)
    100100
    101 !----mode 2 growth 
     101!----mode 2 growth
    102102     typeP = -9;  mode2 = .FALSE.
    103103     do j=1,nz
    … …  
    111111             & zdepthE(k)<z(typeP) .and. &
    112112             & z(typeP)-zdepthE(k)>2*(z(typeP)-z(typeP-1))) then  ! trick that avoids oscillations
    113            deltaz = -avdrhoP(k)/z(typeP)*18./8314.*B  ! conservation of mass 
     113           deltaz = -avdrhoP(k)/z(typeP)*18./8314.*B  ! conservation of mass
    114114           if (deltaz>z(typeP)-z(typeP-1)) then  ! also implies avdrhoP<0.
    115115              mode2 = .TRUE.
    … …  
    143143!
    144144!  Tb<0 initializes temperatures
    145 !  Tb>0 initializes temperatures with Tb 
     145!  Tb>0 initializes temperatures with Tb
    146146!***********************************************************************
    147147  use fast_subs_univ, only: depths_avmeth, equildepth
    … …  
    158158  real(8), intent(INOUT) :: Tb(nz)
    159159  integer, intent(OUT) :: typeF  ! index of depth below which filling occurs
    160   real(8), intent(OUT) :: zdepthE, zdepthF 
     160  real(8), intent(OUT) :: zdepthE, zdepthF
    161161  real(8), intent(OUT) :: ypp(nz) ! (d rho/dt)/Diff
    162162  real(8), intent(OUT) :: Tmean3, zdepthG
    … …  
    171171  real(8) Told(nz), Fsurfold, Tsurfold, Tmean0, avrho2
    172172  real(8) rhosatav0, rhosatav(nz), rlow
    173  
     173
    174174  tmax = EQUILTIME*sols_per_my
    175175  nsteps = int(tmax/dt)     ! calculate total number of timesteps
    176176
    177   Tco2frost = tfrostco2(patm) 
     177  Tco2frost = tfrostco2(patm)
    178178
    179179  !if (Tb<=0.) then  ! initialize
    … …  
    182182     !Tmean0 = Tmean0-5.
    183183     !write(34,*) '# initialized with temperature estimate at',latitude/d2r,'of',Tmean0,'K'
    184      !T(1:nz) = Tmean0 
     184     !T(1:nz) = Tmean0
    185185  !else
    186186     T(1:nz) = Tb
    187187     ! not so good when Fgeotherm is on
    188188  !endif
    189  
     189
    190190  albedo = albedo0
    191191  emiss = emiss0
    … …  
    207207!  Qn=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
    208208  !Qn = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)
    209   !----loop over time steps 
     209  !----loop over time steps
    210210  do n=0,nsteps-1
    211      time = (n+1)*dt         !   time at n+1 
     211     time = (n+1)*dt         !   time at n+1
    212212   !  tdays = time*(sec_per_sol/earthDay) ! parenthesis may improve roundoff
    213213   !  call generalorbit(tdays,a,ecc,omega,eps,marsLs,marsDec,marsR)
    … …  
    215215!     Qnp1=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
    216216     !Qnp1 = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)
    217      
     217
    218218     Tsurfold = Tsurf
    219219     Fsurfold = Fsurf
    … …  
    226226        T(1:nz) = Told(1:nz)
    227227        !call conductionT(nz,z,dt*sec_per_sol,T,Tsurfold,Tco2frost,ti, &
    228              !&              rhocv,Fgeotherm,Fsurf) 
     228             !&              rhocv,Fgeotherm,Fsurf)
    229229        Tsurf = Tco2frost
    230230    !    dE = (- Qn - Qnp1 + Fsurfold + Fsurf + &
    … …  
    240240     endif
    241241     !Qn=Qnp1
    242      
     242
    243243     if (time>=tmax-sols_per_my) then
    244244        Tmean1 = Tmean1 + Tsurf
    … …  
    253253
    254254  enddo  ! end of time loop
    255  
     255
    256256  avrho1 = avrho1/nm
    257257  if (present(avrho1prescribed)) then
    … …  
    266266  endif
    267267  avdrho = avrho2-avrho1
    268   typeP = -9 
     268  typeP = -9
    269269  do i=1,nz
    270270     if (porefill(i)>0.) then
    … …  
    275275  if (typeP>0) then
    276276     avdrhoP = rhosatav(typeP) - avrho1
    277   else 
     277  else
    278278     avdrhoP = -9999.
    279279  end if
    … …  
    282282
    283283  if (Fgeotherm>0.) then
    284      Tb = Tmean1 
     284     Tb = Tmean1
    285285     typeG = 1   ! will be overwritten by depths_avmeth
    286286     rlow = 2*rhosatav(nz)-rhosatav(nz-1)
    … …  
    328328      psv = exp(A/T+B)  ! Clapeyron
    329329
    330 !-----parametrization 3     
     330!-----parametrization 3
    331331!     eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
    332332!     psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)
    333      
     333
    334334      end function psv
    335335

  • trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90

    r3527 r3532  
    3232  real(8) equildepth  ! =zdepthE
    3333  !real(8), external :: zint  ! defined in allinterfaces.mod
    34  
     34
    3535  typeE = -9; equildepth = -9999.
    36   do i=1,nz 
     36  do i=1,nz
    3737     if (rhosatav(i) <= avrho1) then
    3838        typeE=i
    … …  
    6767  integer, intent(OUT) :: typeF  ! index of depth below which filling occurs
    6868  real(8), intent(INOUT) :: zdepthF
    69   real(8), intent(IN) :: B 
     69  real(8), intent(IN) :: B
    7070  real(8), intent(OUT) :: ypp(nz), zdepthG
    7171  integer, intent(INOUT) :: typeG  ! positive on input when Fgeotherm>0
    … …  
    112112
    113113!-depth to shallowest perennial ice
    114   typeP = -9 
     114  typeP = -9
    115115  do i=1,nz
    116116     if (porefill(i)>0.) then
    … …  
    157157  endif
    158158  if (typeG>0 .and. typeT<0) then
    159      call colint(porefill(:)/eta(:),z,nz,typeG-1,nz,cumfillabove) 
     159     call colint(porefill(:)/eta(:),z,nz,typeG-1,nz,cumfillabove)
    160160     newtypeG = -9
    161161     do i=typeG,nz
    … …  
    178178     if (newtypeG>0) typeG=newtypeG
    179179  end if
    180   ! if typeG>0, then all ice at and below typeG should be erased 
     180  ! if typeG>0, then all ice at and below typeG should be erased
    181181end subroutine depths_avmeth
    182182
    … …  
    204204
    205205  ! advance ice table, avdrho>0 is retreat
    206   if (zdepthT>=0. .and. avdrho>0.) then 
     206  if (zdepthT>=0. .and. avdrho>0.) then
    207207     typeP=-9999; typeT=-9999
    208208     do j=1,nz
    … …  
    226226  endif
    227227  if (zdepthT>z(nz)) zdepthT=-9999.
    228  
     228
    229229  ! advance interface, avdrhoP>0 is loss from zdepthP
    230230  if (avdrhoP>0.) then
    … …  
    232232     do j=1,nz
    233233        if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF
    234         if (zdepthT>=0. .and. z(j)>zdepthT) exit 
     234        if (zdepthT>=0. .and. z(j)>zdepthT) exit
    235235        call colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j,integ)
    236236        erase = j
    … …  
    241241
    242242  ! new depth
    243   newtypeP = -9 
     243  newtypeP = -9
    244244  do j=1,nz
    245245     if (zdepthT>=0. .and. z(j)>zdepthT) exit
    … …  
    253253  ub = typeF
    254254  if (newtypeP>0 .and. typeF>0 .and. newtypeP<ub) ub=newtypeP
    255   if (ub>0) then 
     255  if (ub>0) then
    256256     do j=ub,nz
    257257        porefill(j) = porefill(j) + B*ypp(j)

  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3527 r3532  
    501501- Cleaning of "glaciers_mod.F90";
    502502- Optimization for the computation of ice density according to temperature by using a function.
     503
     504== 04/12/2024 == JBC
     505Removing unecessary module/subroutine "interpol_TI_PEM2PCM.F90" + Few small corrections/cleanings.

  • trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90

    r3525 r3532  
    44!-----------------------------------------------------------------------
    55!  Author: LL
    6 !  Purpose: This module gathers the different routines used in the PEM to compute the soil temperature evolution and initialisation. 
    7 ! 
    8 !  Note: depths of layers and mid-layers, soil thermal inertia and 
     6!  Purpose: This module gathers the different routines used in the PEM to compute the soil temperature evolution and initialisation.
     7!
     8!  Note: depths of layers and mid-layers, soil thermal inertia and
    99!        heat capacity are commons in comsoil_PEM.h
    1010!-----------------------------------------------------------------------
    … …  
    1212!=======================================================================
    1313
    14 
    15 
    1614SUBROUTINE compute_tsoil_pem(ngrid,nsoil,firstcall,therm_i,timestep,tsurf,tsoil)
    1715
    … …  
    2422!  Author: LL
    2523!  Purpose: Compute soil temperature using an implict 1st order scheme
    26 ! 
    27 !  Note: depths of layers and mid-layers, soil thermal inertia and 
     24!
     25!  Note: depths of layers and mid-layers, soil thermal inertia and
    2826!        heat capacity are commons in comsoil_PEM.h
    2927!-----------------------------------------------------------------------
    … …  
    3129#include "dimensions.h"
    3230
    33 !-----------------------------------------------------------------------
    34 !  arguments
    35 !  ---------
    36 !  inputs:
    37 integer,                      intent(in) :: ngrid     ! number of (horizontal) grid-points
    38 integer,                      intent(in) :: nsoil     ! number of soil layers 
    39 logical,                      intent(in) :: firstcall ! identifier for initialization call
     31! Inputs:
     32! -------
     33integer,                      intent(in) :: ngrid     ! number of (horizontal) grid-points
     34integer,                      intent(in) :: nsoil     ! number of soil layers
     35logical,                      intent(in) :: firstcall ! identifier for initialization call
    4036real, dimension(ngrid,nsoil), intent(in) :: therm_i   ! thermal inertia [SI]
    4137real,                         intent(in) :: timestep  ! time step [s]
    4238real, dimension(ngrid),       intent(in) :: tsurf     ! surface temperature [K]
    43  
    44 ! outputs:
     39! Outputs:
     40!---------
    4541real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
    46 ! local variables:
    47 integer :: ig, ik   
     42! Local:
     43!-------
     44integer :: ig, ik
    4845
    4946! 0. Initialisations and preprocessing step
    … …  
    5249    do ig = 1,ngrid
    5350        do ik = 0,nsoil - 1
    54             mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa   
     51            mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa
    5552        enddo
    5653    enddo
    … …  
    6461    enddo
    6562
    66 ! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
     63! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEM
    6764    ! mu_PEM
    6865    mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
    … …  
    9996! Other layers:
    10097        do ik = 1,nsoil - 1
    101                 tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik) 
     98                tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
    10299        enddo
    103100    enddo
    … …  
    139136#include "dimensions.h"
    140137
    141 !-----------------------------------------------------------------------
    142 !  arguments
    143 !  ---------
    144 !  inputs:
     138! Inputs:
     139!--------
    145140integer,                      intent(in) :: ngrid   ! number of (horizontal) grid-points
    146141integer,                      intent(in) :: nsoil   ! number of soil layers
    147142real, dimension(ngrid,nsoil), intent(in) :: therm_i ! thermal inertia [SI]
    148143real, dimension(ngrid),       intent(in) :: tsurf   ! surface temperature [K]
    149 
    150 ! outputs:
     144! Outputs:
     145!---------
    151146real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
    152 ! local variables:
     147! Local:
     148!-------
    153149integer :: ig, ik, iloop
    154150
    … …  
    169165enddo
    170166
    171 ! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
     167! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEM
    172168! mu_PEM
    173169mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
    … …  
    179175do ig = 1,ngrid
    180176    ! d_k
    181     do ik = 1,nsoil-1
     177    do ik = 1,nsoil - 1
    182178        coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
    183179    enddo
    … …  
    207203!  2. Compute soil temperatures
    208204do iloop = 1,10 !just convergence
    209     ! First layer:
    210     do ig = 1,ngrid
    211         tsoil(ig,1)=(tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
    212                    (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
    213     ! Other layers:
    214     do ik = 1,nsoil - 1
    215         tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
    216     enddo
    217 enddo
    218 
     205    do ig = 1,ngrid
     206        ! First layer:
     207        tsoil(ig,1) = (tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
     208                      (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
     209        ! Other layers:
     210        do ik = 1,nsoil - 1
     211            tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
     212        enddo
     213    enddo
    219214enddo ! iloop
    220215

  • trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90

    r3161 r3532  
    33implicit none
    44
    5 integer, parameter                        :: nsoilmx_PEM = 68   ! number of layers in the PEM
    6 real, allocatable, dimension(:),     save :: layer_PEM          ! soil layer depths [m]
    7 real, allocatable, dimension(:),     save :: mlayer_PEM         ! soil mid-layer depths [m]
    8 real, allocatable, dimension(:,:,:), save :: TI_PEM             ! soil thermal inertia [SI]
    9 real, allocatable, dimension(:,:),   save :: inertiedat_PEM     ! soil thermal inertia saved as reference for current climate [SI]
     5integer, parameter                  :: nsoilmx_PEM = 68   ! number of layers in the PEM
     6real, allocatable, dimension(:)    :: layer_PEM          ! soil layer depths [m]
     7real, allocatable, dimension(:)    :: mlayer_PEM         ! soil mid-layer depths [m]
     8real, allocatable, dimension(:,:,:) :: TI_PEM             ! soil thermal inertia [SI]
     9real, allocatable, dimension(:,:)  :: inertiedat_PEM     ! soil thermal inertia saved as reference for current climate [SI]
    1010! Variables (FC: built in firstcall in soil.F)
    11 real, allocatable, dimension(:,:,:), save :: tsoil_PEM          ! sub-surface temperatures [K]
    12 real, allocatable, dimension(:,:),   save :: mthermdiff_PEM     ! (FC) mid-layer thermal diffusivity [SI]
    13 real, allocatable, dimension(:,:),   save :: thermdiff_PEM      ! (FC) inter-layer thermal diffusivity [SI]
    14 real, allocatable, dimension(:),     save :: coefq_PEM          ! (FC) q_{k+1/2} coefficients [SI]
    15 real, allocatable, dimension(:,:),   save :: coefd_PEM          ! (FC) d_k coefficients [SI]
    16 real, allocatable, dimension(:,:),   save :: alph_PEM           ! (FC) alpha_k coefficients [SI]
    17 real, allocatable, dimension(:,:),   save :: beta_PEM           ! beta_k coefficients [SI]
    18 real,                                save :: mu_PEM             ! mu coefficient [SI]
    19 real,                                save :: fluxgeo            ! Geothermal flux [W/m^2]
    20 real,                                save :: depth_breccia      ! Depth at which we have breccia [m]
    21 real,                                save :: depth_bedrock      ! Depth at which we have bedrock [m]
    22 integer,                             save :: index_breccia      ! last index of the depth grid before having breccia
    23 integer,                             save :: index_bedrock      ! last index of the depth grid before having bedrock
    24 logical                                   :: soil_pem           ! True by default, to run with the subsurface physic. Read in pem.def
    25 real,                                save :: ini_huge_h2oice    ! Initial value for huge reservoirs of H2O ice [kg/m^2]
    26 logical,                             save :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure
     11real, allocatable, dimension(:,:,:) :: tsoil_PEM          ! sub-surface temperatures [K]
     12real, allocatable, dimension(:,:)  :: mthermdiff_PEM     ! (FC) mid-layer thermal diffusivity [SI]
     13real, allocatable, dimension(:,:)  :: thermdiff_PEM      ! (FC) inter-layer thermal diffusivity [SI]
     14real, allocatable, dimension(:)    :: coefq_PEM          ! (FC) q_{k+1/2} coefficients [SI]
     15real, allocatable, dimension(:,:)  :: coefd_PEM          ! (FC) d_k coefficients [SI]
     16real, allocatable, dimension(:,:)  :: alph_PEM           ! (FC) alpha_k coefficients [SI]
     17real, allocatable, dimension(:,:)  :: beta_PEM           ! beta_k coefficients [SI]
     18real                                :: mu_PEM             ! mu coefficient [SI]
     19real                                :: fluxgeo            ! Geothermal flux [W/m^2]
     20real                                :: depth_breccia      ! Depth at which we have breccia [m]
     21real                                :: depth_bedrock      ! Depth at which we have bedrock [m]
     22integer                            :: index_breccia      ! last index of the depth grid before having breccia
     23integer                            :: index_bedrock      ! last index of the depth grid before having bedrock
     24logical                             :: soil_pem           ! True by default, to run with the subsurface physic. Read in pem.def
     25real                                :: ini_huge_h2oice    ! Initial value for huge reservoirs of H2O ice [kg/m^2]
     26logical                            :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure
    2727
    2828!=======================================================================

  • trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90

    r3527 r3532  
    2424!!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do
    2525!!!          the ice transfer
    26 !!!         
    27 !!!         
     26!!!
     27!!!
    2828!!! Author: LL
    2929!!!
    … …  
    4747real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
    4848real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
    49 ! Local 
     49! Local
    5050!------
    5151real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
    52 real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow 
     52real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow
    5353
    5454write(*,*) "Flow of CO2 glaciers"
    … …  
    6767!!! Purpose: Main for H2O glaciers evolution: compute maximum thickness, and do
    6868!!!          the ice transfer
    69 !!!         
    70 !!!         
     69!!!
     70!!!
    7171!!! Author: LL
    7272!!!
    … …  
    8787real, dimension(ngrid,nslope), intent(inout) :: flag_h2oflow      ! flag to see if there is flow on the subgrid slopes
    8888real, dimension(ngrid),        intent(inout) :: flag_h2oflow_mesh ! same but within the mesh
    89 ! Local 
    90 real, dimension(ngrid,nslope) :: hmax ! Physical x Slope field: maximum thickness for co2  glacier before flow 
     89! Local
     90real, dimension(ngrid,nslope) :: hmax ! Physical x Slope field: maximum thickness for co2  glacier before flow
    9191
    9292write(*,*) "Flow of H2O glaciers"
    … …  
    111111!!!
    112112!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle before initating flow
    113 !!!         
     113!!!
    114114!!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
    115115!!!
    … …  
    222222            endif ! co2ice > hmax
    223223        endif ! iflat
    224     enddo !islope 
     224    enddo !islope
    225225enddo !ig
    226226
    … …  
    237237!!!
    238238!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    239  
     239
    240240use constants_marspem_mod, only: alpha_clap_co2,beta_clap_co2
    241241
    242 implicit none 
     242implicit none
    243243
    244244! arguments:
    … …  
    253253
    254254! OUTPUT
    255 real, dimension(ngrid,nslope), intent(out) :: Tcond ! Physical points: condensation temperature of CO2, yearly averaged 
     255real, dimension(ngrid,nslope), intent(out) :: Tcond ! Physical points: condensation temperature of CO2, yearly averaged
    256256
    257257! LOCAL

  • trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90

    r3527 r3532  
    401401! Locals
    402402! ------
    403 integer :: ik       ! Loop variables 
     403integer :: ik       ! Loop variables
    404404integer :: indexice ! Index of the ice
    405405
    … …  
    412412    if(ice_depth < mlayer_PEM(0)) then
    413413        indexice = 0.
    414     else     
     414    else
    415415        do ik = 0,nsoil - 2 ! go through all the layers to find the ice locations
    416416            if (mlayer_PEM(ik) <= ice_depth .and. mlayer_PEM(ik + 1) > ice_depth) then

  • trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90

    r3512 r3532  
    3131logical       :: ok
    3232integer       :: cstat
    33 character(10) :: frmt
     33character(50) :: info_frmt
    3434character(20) :: ich1, ich2, ich3, ich4, fch1, fch2, fch3
    3535
    … …  
    5353    ! Martian date, Number of Martians years done by the PEM run, Number of Martians years done by the chainded simulation, Code of the stopping criterion
    5454    ! The conversion ratio from Planetary years to Earth years is given in the header of the file
    55     write(1,*) year_bp_ini + i_myear, i_myear_leg, i_myear, stopPEM
     55    info_frmt = '(f'//int2str(nb_digits(year_bp_ini + i_myear) + 5)//'.4,'//'f'//int2str(nb_digits(i_myear_leg) + 5)//'.4,'//'f'//int2str(nb_digits(i_myear) + 5)//'.4,'//'i0)'
     56    write(1,trim(adjustl(info_frmt))) year_bp_ini + i_myear, i_myear_leg, i_myear, stopPEM
    5657    close(1)
    5758else

  • trunk/LMDZ.COMMON/libf/evolution/math_mod.F90

    r3526 r3532  
    33!!!
    44!!! Purpose: The module contains all the mathematical SUBROUTINE used in the PEM
    5 !!!           
     5!!!
    66!!! Author: Adapted from Schorgofer MSIM (N.S, Icarus 2010), impletented here by LL
    77!!!
    … …  
    1717!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1818!!!
    19 !!! Purpose: Compute the first derivative of a function y(z) on an irregular grid 
    20 !!!           
     19!!! Purpose: Compute the first derivative of a function y(z) on an irregular grid
     20!!!
    2121!!! Author: From N.S (N.S, Icarus 2010), impletented here by LL
    2222!!!
    … …  
    5959!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    6060!!!
    61 !!! Purpose: Compute the second derivative of a function y(z) on an irregular grid 
    62 !!!           
     61!!! Purpose: Compute the second derivative of a function y(z) on an irregular grid
     62!!!
    6363!!! Author: N.S (raw copy/paste from MSIM)
    6464!!!
    … …  
    9797!!!
    9898!!! Purpose: First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
    99 !!!           
     99!!!
    100100!!! Author: N.S (raw copy/paste from MSIM)
    101101!!!
    … …  
    128128!!!
    129129!!! Purpose:  Column integrates y on irregular grid
    130 !!!           
     130!!!
    131131!!! Author: N.S (raw copy/paste from MSIM)
    132132!!!
    133133!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    134  
     134
    135135implicit none
    136  
     136
    137137integer, intent(in) :: nz, i1, i2
    138138real,    intent(in) :: y(nz), z(nz)
    … …  
    156156!!!
    157157!!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
    158 !!! 
     158!!!
    159159!!! Author: LL
    160160!!!

  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r3527 r3532  
    6262use compute_tend_mod,           only: compute_tend
    6363use info_PEM_mod,               only: info_PEM
    64 use interpol_TI_PEM2PCM_mod,    only: interpol_TI_PEM2PCM
    6564use nb_time_step_PCM_mod,       only: nb_time_step_PCM
    6665use pemetat0_mod,               only: pemetat0
    … …  
    298297    call execute_command_line('squeue -j '//trim(job_id)//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat)
    299298    if (cstat /= 0) then
    300         call execute_command_line('qstat -f '//trim(job_id)//' | grep "Walltime" | awk ''{print $3}'' > tmp_cmdout.txt', cmdstat = cstat)
     299        call execute_command_line('qstat -f '//trim(job_id)//' | grep "Walltime" | awk ''{print $3}'' > tmp_cmdout.txt',cmdstat = cstat)
    301300        if (cstat > 0) then
    302301            error stop 'pem: command execution failed!'
    … …  
    939938            do ig = 1,ngrid
    940939                do islope = 1,nslope
    941                     call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps(ig),sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2),ice_porefilling(ig,:,islope),porefill,ssi_depth)
     940                    call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps(ig),(/sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2)/),ice_porefilling(ig,:,islope),porefill,ssi_depth)
    942941                    icetable_depth(ig,islope) = ssi_depth
    943942                    ice_porefilling(ig,:,islope) = porefill
    … …  
    961960            do ig = 1,ngrid
    962961                do islope = 1,nslope
    963                     do l = 1,nsoilmx_PEM 
     962                    do l = 1,nsoilmx_PEM
    964963                        if (l == 1) then
    965964                            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
    … …  
    10081007                call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
    10091008                call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope))
    1010             endif                       
     1009            endif
    10111010        endif
    10121011    enddo
    … …  
    11011100! III_a.2 Tsoil update (for startfi)
    11021101if (soil_pem) then
    1103     call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
     1102    inertiesoil = TI_PEM(:,:nsoilmx,:)
    11041103    tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom(:,:,:,timelen)
    11051104#ifndef CPP_STD

  • trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90

    r3525 r3532  
    4545real,                           intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa]
    4646real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr1       ! surface temperature at the first year of PCM call [K]
    47 real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2       ! surface temperature at the second  year of PCM call [K]
     47real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2       ! surface temperature at the second year of PCM call [K]
    4848real, dimension(ngrid,timelen), intent(in) :: q_co2               ! MMR tracer co2 [kg/kg]
    4949real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
    … …  
    281281                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
    282282            else
    283 ! predictor corrector: restart from year 1 of the PCM and build the evolution of
    284 ! tsoil at depth
     283! predictor corrector: restart from year 1 of the PCM and build the evolution of tsoil at depth
    285284                tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
    286285                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
    … …  
    289288                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
    290289
    291                 do iloop = nsoil_PCM+1,nsoil_PEM
     290                do iloop = nsoil_PCM + 1,nsoil_PEM
    292291                    tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
    293292                enddo
    … …  
    295294
    296295            do it = 1,timelen
    297                 do isoil = nsoil_PCM+1,nsoil_PEM
    298                     tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
    299                 enddo
    300             enddo
    301             do isoil = nsoil_PCM+1,nsoil_PEM
    302                 do ig = 1,ngrid
    303                     watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
    304                 enddo
     296                tsoil_inst(:,nsoil_PCM + 1:nsoil_PEM,islope,it) = tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)
     297            enddo
     298            do isoil = nsoil_PCM + 1,nsoil_PEM
     299                watersoil_avg(:,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(:,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(:,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
    305300            enddo
    306301        enddo ! islope
    … …  
    314309                write(*,*)'PEM settings: failed loading <ice_table_depth>'
    315310                write(*,*)'will reconstruct the values of the ice table given the current state'
    316                 call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
     311                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
    317312                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
    318313                do islope = 1,nslope
    … …  
    322317            write(*,*) 'PEMETAT0: ICE TABLE done'
    323318        else if (icetable_dynamic) then
     319            call get_field("ice_porefilling",ice_porefilling,found)
     320            if (.not. found) then
     321                write(*,*)'PEM settings: failed loading <ice_porefilling>'
     322                ice_porefilling = 0.
     323            endif
    324324            call get_field("ice_table_depth",ice_table_depth,found)
    325325            if (.not. found) then
    … …  
    331331                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
    332332                enddo
    333             endif
    334             call get_field("ice_porefilling",ice_porefilling,found)
    335             if (.not. found) then
    336                 write(*,*)'PEM settings: failed loading <ice_porefilling>'
    337                 ice_porefilling = 1.
    338333            endif
    339334            write(*,*) 'PEMETAT0: ICE TABLE done'
    … …  
    510505            write(*,*) 'PEMETAT0: Ice table done'
    511506        else if (icetable_dynamic) then
    512             ice_porefilling = 1.
     507            ice_porefilling = 0.
    513508            ice_table_depth = -9999.
    514509            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
    … …  
    518513            write(*,*) 'PEMETAT0: Ice table done'
    519514        endif
    520        
     515
    521516!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    522517!d) Regolith adsorbed

  • trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90

    r3498 r3532  
    4141
    4242write(*,*) "Update of the CO2 tendency from the current pressure"
    43    
     43
    4444! Evolution of the water ice for each physical point
    4545do i = 1,ngrid

  • trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90

    r3206 r3532  
    1515
    1616!=======================================================================
    17 !  Author: LL, based on work by  Ehouarn Millour (07/2006)
     17!  Author: LL, based on work by Ehouarn Millour (07/2006)
    1818!
    1919!  Purpose: Read and/or initialise soil depths and properties
    … …  
    5353alpha = 2
    5454do iloop = 0,nsoil_PCM - 1
    55     mlayer_PEM(iloop) = lay1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
     55    mlayer_PEM(iloop) = lay1*(1 + iloop**2.9*(1 - exp(-real(iloop)/20.)))
    5656enddo
    5757
    58 do iloop = 0, nsoil_PEM-1
    59     if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM-1)) then
     58do iloop = 0,nsoil_PEM - 1
     59    if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM - 1)) then
    6060        index_powerlaw = iloop
    6161        exit
    … …  
    6363enddo
    6464
    65 do iloop = nsoil_PCM, nsoil_PEM-1
    66     mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop-nsoil_PCM)-0.5))
     65do iloop = nsoil_PCM,nsoil_PEM - 1
     66    mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop - nsoil_PCM) - 0.5))
    6767enddo
    6868
    6969! Build layer_PEM()
    70 do iloop=1,nsoil_PEM-1
    71 layer_PEM(iloop)=(mlayer_PEM(iloop)+mlayer_PEM(iloop-1))/2.
     70do iloop = 1,nsoil_PEM - 1
     71layer_PEM(iloop) = (mlayer_PEM(iloop) + mlayer_PEM(iloop - 1))/2.
    7272enddo
    73 layer_PEM(nsoil_PEM)=2*mlayer_PEM(nsoil_PEM-1) - mlayer_PEM(nsoil_PEM-2)
     73layer_PEM(nsoil_PEM) = 2*mlayer_PEM(nsoil_PEM - 1) - mlayer_PEM(nsoil_PEM - 2)
    7474
    7575

  • trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90

    r3512 r3532  
    680680    call abort_physic("writediagsoilpem","firstname too short",1)
    681681  endif
    682  
     682
    683683  ! Set output sample rate
    684684  isample = int(ecritpem) ! same as for diagpem outputs
    685  
     685
    686686  ! Create output NetCDF file
    687687  if (is_master) then
    … …  
    724724    enddo
    725725   endif
    726    
     726
    727727   ! write "header" of file (longitudes, latitudes, geopotential, ...)
    728728   if (klon_glo>1) then ! general 3D case
    … …  
    733733
    734734  endif ! of if (is_master)
    735  
     735
    736736  ! set zitau to -1 to be compatible with zitau incrementation step below
    737737  zitau=-1
    738  
     738
    739739else
    740740  ! If not an initialization call, simply open the NetCDF file
    … …  
    757757    ntime=ntime+1
    758758    date = float(zitau + 1)
    759    
     759
    760760    if (is_master) then
    761761     ! Get NetCDF ID for "time"
    … …  
    769769     if (ierr.ne.NF_NOERR) then
    770770      write(*,*)"writediagsoilpem: Failed writing date to time variable"
    771       call abort_physic("writediagsoilpem","failed writing time",1) 
     771      call abort_physic("writediagsoilpem","failed writing time",1)
    772772     endif
    773773    endif ! of if (is_master)
    … …  
    810810   data3_1d(1,1:nsoilmx_PEM)=px(1,1:nsoilmx_PEM)
    811811  endif
    812  
     812
    813813  ! B. Write (append) the variable to the NetCDF file
    814814  if (is_master) then
    … …  
    828828    call def_var(nid,name,title,units,4,id,varid,ierr)
    829829  endif ! of if (ierr.ne.NF_NOERR)
    830  
     830
    831831  ! B.2. Prepare things to be able to write/append the variable
    832832  corners(1)=1
    … …  
    834834  corners(3)=1
    835835  corners(4)=ntime
    836  
     836
    837837  if (klon_glo==1) then
    838838    edges(1)=1
    … …  
    843843  edges(3)=nsoilmx_PEM
    844844  edges(4)=1
    845  
     845
    846846  ! B.3. Write the slab of data
    847847!#ifdef NC_DOUBLE
    … …  
    917917  corners(2)=1
    918918  corners(3)=ntime
    919  
     919
    920920  if (klon_glo==1) then
    921921    edges(1)=1
    … …  
    925925  edges(2)=nbp_lat
    926926  edges(3)=1
    927  
     927
    928928  ! B.3. Write the slab of data
    929929!#ifdef NC_DOUBLE
    … …  
    966966  ! B.2. Prepare things to be able to write/append the variable
    967967  corners(1)=ntime
    968  
     968
    969969  edges(1)=1
    970970
Note: See TracChangeset for help on using the changeset viewer.