Changeset 3526 for trunk


Ignore:
Timestamp:
Nov 19, 2024, 5:54:18 PM (2 days ago)
Author:
jbclement
Message:

PEM:
Removing unused or redundant Norbert Schorghofer's subroutines (follow-up of r3493) + cleaning and some modifications of related subroutines.
JBC

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

Legend:

Unmodified
Added
Removed

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

    r3512 r3526  
    1 
    21!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    32!!!
    … …  
    171170 
    172171end subroutine dyn_ss_ice_m
    173 
    174 
    175 

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

    r3493 r3526  
    1 
    21module allinterfaces
    32  ! interfaces from Fortran 90 subroutines and functions
    … …  
    8685       real(8), intent(OUT) :: ypp(nz), zdepthG
    8786       integer, intent(INOUT) :: typeG
    88        real(8), external :: zint, deriv1_onesided, colint
     87       real(8), external :: zint
    8988     end subroutine depths_avmeth
    9089  end interface
    … …  
    9695       real(8) constriction
    9796     end function constriction
    98   end interface
    99  
    100   interface
    101      subroutine assignthermalproperties(nz,thIn,rhoc,ti,rhocv, &
    102           &                      typeT,icefrac,porosity,porefill)
    103        implicit none
    104        integer, intent(IN) :: nz
    105        integer, intent(IN), optional :: typeT
    106        real(8), intent(IN), optional :: icefrac
    107        real(8), intent(IN) :: thIn, rhoc
    108        real(8), intent(IN), optional :: porosity, porefill(nz)
    109        real(8), intent(OUT) :: ti(nz), rhocv(nz)
    110      end subroutine assignthermalproperties
    111   end interface
    112 
    113   interface
    114      pure subroutine icechanges_poreonly(nz,z,typeF,typeG,avdrhoP,ypp,B,porefill)
    115        implicit none
    116        integer, intent(IN) :: nz, typeF, typeG
    117        real(8), intent(IN) :: z(nz), ypp(nz), avdrhoP, B
    118        real(8), intent(INOUT) :: porefill(nz)
    119      end subroutine icechanges_poreonly
    12097  end interface
    12198
    … …  
    131108  end interface
    132109
    133   interface
    134      subroutine compactoutput(unit,porefill,nz)
    135        implicit none
    136        integer, intent(IN) :: unit,nz
    137        real(8), intent(IN) :: porefill(nz)
    138      end subroutine compactoutput
    139   end interface
    140 
    141110  !end of fast_subs_univ
    142   !begin derivs.f90
    143 
    144   interface
    145      subroutine deriv1(z,nz,y,y0,yNp1,yp)
    146        implicit none
    147        integer, intent(IN) :: nz
    148        real(8), intent(IN) :: z(nz),y(nz),y0,yNp1
    149        real(8), intent(OUT) :: yp(nz)
    150      end subroutine deriv1
    151   end interface
    152 
    153   interface
    154      subroutine deriv2_full(z,nz,a,b,a0,b0,bNp1,yp2)
    155        implicit none
    156        integer, intent(IN) :: nz
    157        real(8), intent(IN) :: z(nz),a(nz),b(nz),a0,b0,bNp1
    158        real(8), intent(OUT) :: yp2(nz)
    159      end subroutine deriv2_full
    160   end interface
    161 
    162   interface
    163      subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2)
    164        implicit none
    165        integer, intent(IN) :: nz
    166        real(8), intent(IN) :: z(nz),y(nz),y0,yNp1
    167        real(8), intent(OUT) :: yp2(nz)
    168      end subroutine deriv2_simple
    169   end interface
    170   interface
    171      real(8) pure function deriv1_onesided(j,z,nz,y)
    172        implicit none
    173        integer, intent(IN) :: nz,j
    174        real(8), intent(IN) :: z(nz),y(nz)
    175      end function deriv1_onesided
    176   end interface
    177 
    178   !end of derivs.f90
    179   !begin fast_subs_exper.f90
    180 
    181   interface
    182      subroutine icelayer_exper(bigstep, nz, thIn, rhoc, z, porosity, pfrost, &
    183           & zdepthT, icefrac, zdepthF, zdepthE, porefill, Tmean, Tampl, Diff, zdepthG)
    184        use miscparameters, only : d2r, NMAX, icedensity
    185        implicit none
    186        integer, intent(IN) :: nz
    187        real(8), intent(IN) :: bigstep
    188        real(8), intent(IN) :: thIn, rhoc, z(NMAX), porosity, pfrost
    189        real(8), intent(INOUT) :: zdepthT, zdepthF, porefill(nz)
    190        real(8), intent(OUT) :: zdepthE
    191        real(8), intent(IN) :: icefrac, Diff, Tmean, Tampl
    192        real(8), intent(OUT) :: zdepthG
    193      end subroutine icelayer_exper
    194   end interface
    195 
    196   interface
    197      subroutine ajsub_exper(typeT, nz, z, ti, rhocv, pfrost, Tmean, Tampl, &
    198           &     avdrho, avdrhoP, zdepthE, typeF, zdepthF, ypp, porefill, &
    199           &     B, typeG, zdepthG)
    200        use miscparameters, only : NMAX, solsperyear, marsDay
    201        implicit none
    202        integer, intent(IN) :: nz, typeT
    203        real(8), intent(IN) :: z(NMAX), ti(NMAX), rhocv(NMAX), pfrost
    204        real(8), intent(IN) :: Tmean, Tampl
    205        real(8), intent(OUT) :: avdrho, avdrhoP
    206        real(8), intent(OUT) :: zdepthE
    207        integer, intent(OUT) :: typeF
    208        real(8), intent(INOUT) :: zdepthF
    209        real(8), intent(OUT) :: ypp(nz)
    210        real(8), intent(IN) :: porefill(nz)
    211        real(8), intent(IN) :: B
    212        integer, intent(OUT) :: typeG
    213        real(8), intent(OUT) :: zdepthG
    214        real(8), external :: Tsurface, psv
    215        real(8), external :: equildepth
    216      end subroutine ajsub_exper
    217   end interface
    218  
    219   !end fast_subs_exper
    220 
    221111  ! Other
    222   interface
    223      pure function flux_mars77(R,decl,HA,albedo,fracir,fracscat)
    224        implicit none
    225        real*8 flux_mars77
    226        real*8, intent(IN) :: R,decl,HA,albedo,fracIR,fracScat
    227      end function flux_mars77
    228   end interface
    229 
    230   interface
    231      pure function colint(y,z,nz,i1,i2)
    232        implicit none
    233        integer, intent(IN) :: nz, i1, i2
    234        real(8), intent(IN) :: y(nz),z(nz)
    235        real(8) colint
    236      end function colint
    237   end interface
    238 
    239112  interface
    240113   subroutine dyn_ss_ice_m(ssi_depth_in,T1,T_in,nsoil,thIn,p0,pfrost,porefill_in,porefill,ssi_depth)
    … …  
    250123  end interface
    251124
    252     interface
    253       subroutine dyn_ss_ice_m_wrapper(ngrid,nsoil,tHIn,p0,pfrost,T_in,ssi_depth_in,porefill_in,porefill,ssi_depth)
    254            implicit none
    255            integer, intent(IN) :: nsoil,ngrid
    256            real(8),  intent(IN) :: thIn(ngrid),ssi_depth_in(ngrid)
    257            real(8),  intent(IN) :: p0(ngrid), pfrost(ngrid)
    258            real(8),  intent(IN) :: T_in(nsoil,ngrid)
    259            real(8), intent(OUT) :: porefill(nsoil,ngrid)
    260            real(8), intent(IN) :: porefill_in(nsoil,ngrid)
    261            real(8), intent(OUT) :: ssi_depth(ngrid)
    262    end subroutine dyn_ss_ice_m_wrapper
    263   end interface
    264 
    265125end module allinterfaces

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

    r3512 r3526  
    294294
    295295
    296 subroutine outputmoduleparameters
    297   use miscparameters
    298   use thermalmodelparam_mars
    299   implicit none
    300 !  print *,'Parameters stored in modules'
    301 !  print *,'  Ice bulk density',icedensity,'kg/m^3'
    302 !  print *,'  dt=',dt,'Mars solar days'
    303 !  print *,'  Fgeotherm=',Fgeotherm,'W/m^2'
    304 !  write(6,'(2x,a27,1x,f5.3)') 'Emissivity of dry surface=',emiss0
    305 !  write(6,'(2x,a12,1x,f5.3,2x,a16,1x,f5.3)') 'CO2 albedo=',co2albedo,'CO2 emissivity=',co2emiss
    306 !  print *,'  Thermal model equilibration time',EQUILTIME,'Mars years'
    307 end subroutine outputmoduleparameters
    308 
    309 
    310 
     296real*8 function tfrostco2(p)
     297!     the inverse of function psvco2
     298!     input is pressure in Pascal, output is temperature in Kelvin
     299implicit none
     300real*8 p
     301tfrostco2 = 3182.48/(23.3494+log(100./p))
     302end function

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

    r3493 r3526  
    5555!***********************************************************************
    5656  use allinterfaces, except_this_one => depths_avmeth
     57  use math_mod, only: deriv2_simple, deriv1_onesided, deriv1, colint
    5758  implicit none
    5859  integer, intent(IN) :: nz, typeT
    … …  
    118119  call deriv1(z,nz,eta(:),1.d0,eta(nz-1),ap)
    119120  if (typeP>0 .and. typeP<nz-2) then
    120      ap_one=deriv1_onesided(typeP,z,nz,eta(:))
     121     call deriv1_onesided(typeP,z,nz,eta(:),ap_one)
    121122     ! print *,typeP,ap(typeP),ap_one
    122123     ap(typeP)=ap_one
    … …  
    151152  endif
    152153  if (typeG>0 .and. typeT<0) then
    153      cumfillabove = colint(porefill(:)/eta(:),z,nz,typeG-1,nz)
     154     call colint(porefill(:)/eta(:),z,nz,typeG-1,nz,cumfillabove)
    154155     newtypeG = -9
    155156     do i=typeG,nz
    156157        if (minval(eta(i:nz))<=0.) cycle  ! eta=0 means completely full
    157         cumfill=colint(porefill(:)/eta(:),z,nz,i,nz)
     158        call colint(porefill(:)/eta(:),z,nz,i,nz,cumfill)
    158159        if (cumfill<yp(i)*18./8314.*B) then  ! usually executes on i=typeG
    159160           if (i>typeG) then
    … …  
    194195
    195196
    196 pure subroutine icechanges_poreonly(nz,z,typeF,typeG,avdrhoP,ypp,B,porefill)
    197   use allinterfaces, except_this_one => icechanges_poreonly
    198   implicit none
    199   integer, intent(IN) :: nz, typeF, typeG
    200   real(8), intent(IN) :: z(nz), ypp(nz), avdrhoP, B
    201   real(8), intent(INOUT) :: porefill(nz)
    202   integer j, erase, newtypeP, ub
    203   real(8) integ
    204  
    205   !----retreat
    206   ! avdrhoP>0 is outward loss from zdepthP
    207   ! avdrhoP<0 means gain at zdepthP or no ice anywhere
    208   if (avdrhoP>0.) then
    209      erase=0
    210      do j=1,nz
    211         if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF
    212         integ = colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j)
    213         erase = j
    214         if (integ>B*avdrhoP*18./8314.) exit
    215      end do
    216      if (erase>0) porefill(1:erase)=0.
    217   endif
    218 
    219   ! new depth
    220   newtypeP = -9
    221   do j=1,nz
    222      if (porefill(j)>0.) then
    223         newtypeP = j  ! first point with ice
    224         exit
    225      endif
    226   enddo
    227 
    228   !----diffusive filling
    229   ub = typeF
    230   if (newtypeP>0 .and. typeF>0 .and. newtypeP<ub) ub=newtypeP
    231   if (ub>0) then 
    232      do j=ub,nz
    233         ! B=Diff/(porosity*icedensity)*86400*365.24*bigstep
    234         porefill(j) = porefill(j) + B*ypp(j)
    235         if (porefill(j)<0.) porefill(j)=0.
    236         if (porefill(j)>1.) porefill(j)=1.
    237      enddo
    238   end if
    239  
    240   !----enact bottom boundary
    241   if (typeG>0) porefill(typeG:nz)=0.
    242  
    243 end subroutine icechanges_poreonly
    244 
    245 
    246 
    247197pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, &
    248198     & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG)
    … …  
    254204  use miscparameters, only : icedensity
    255205  use allinterfaces, except_this_one => icechanges
     206  use math_mod, only: colint
    256207  implicit none
    257208  integer, intent(IN) :: nz, typeF, typeG
    … …  
    293244        if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF
    294245        if (zdepthT>=0. .and. z(j)>zdepthT) exit
    295         integ = colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j)
     246        call colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j,integ)
    296247        erase = j
    297248        if (integ>B*avdrhoP*18./8314.) exit
    … …  
    332283  end if
    333284end subroutine icechanges
    334 
    335 
    336 subroutine assignthermalproperties(nz,thIn,rhoc, &
    337      &    ti,rhocv,typeT,icefrac,porosity,porefill)
    338 !*********************************************************
    339 ! assign thermal properties of soil
    340 !*********************************************************
    341   implicit none
    342   integer, intent(IN) :: nz
    343   integer, intent(IN), optional :: typeT
    344   real(8), intent(IN), optional :: icefrac
    345   real(8), intent(IN) :: thIn, rhoc
    346   real(8), intent(IN), optional :: porosity, porefill(nz)
    347   real(8), intent(OUT) :: ti(nz), rhocv(nz)
    348   integer j
    349   real(8) newrhoc, newti, fill
    350   real(8), parameter :: NULL=0.
    351 
    352   ti(1:nz) = thIn
    353   rhocv(1:nz) = rhoc
    354   if (typeT>0) then
    355      call soilthprop(porosity,NULL,rhoc,thIn,2,newrhoc,newti,icefrac)
    356      rhocv(typeT:nz) = newrhoc
    357      ti(typeT:nz) = newti
    358   endif
    359   do j=1,nz
    360      fill = porefill(j)   ! off by half point
    361      if (fill>0. .and. (typeT<0 .or. (typeT>0 .and. j<typeT))) then
    362         call soilthprop(porosity,fill,rhoc,thIn,1,rhocv(j),ti(j),NULL)
    363      endif
    364   enddo
    365 end subroutine assignthermalproperties
    366 
    367 
    368 
    369 subroutine compactoutput(unit,porefill,nz)
    370   implicit none
    371   integer, intent(IN) :: unit,nz
    372   real(8), intent(IN) :: porefill(nz)
    373   integer j
    374   do j=1,nz
    375      if (porefill(j)==0.) then
    376         write(unit,'(1x,f2.0)',advance='no') porefill(j)
    377      else
    378         write(unit,'(1x,f7.5)',advance='no') porefill(j)
    379      endif
    380   enddo
    381   write(unit,"('')")
    382 end subroutine compactoutput
    383 

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

    r3493 r3526  
    1212  !   length of Earth year in days = 365.24
    1313end module miscparameters
    14 

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

    r3493 r3526  
    7676 
    7777end subroutine soilthprop
    78 
    79 
    80      
    81 subroutine smartgrid(nz,z,zdepth,thIn,rhoc,porosity,ti,rhocv,layertype,icefrac)
    82 !***********************************************************************
    83 ! smartgrid: returns intelligently spaced grid and appropriate
    84 !            values of thermal inertia ti and rho*c in icy layer
    85 !                 
    86 !     INPUTS:
    87 !             nz = number of grid points
    88 !             z = grid spacing for dry regolith
    89 !                 (will be partially overwritten)
    90 !             zdepth = depth where ice table starts
    91 !                      negative values indicate no ice
    92 !             rhoc = heat capacity per volume of dry regolith [J/m^3]
    93 !             thIn = thermal inertia of dry regolith [SI-units]
    94 !             porosity = void space / total volume
    95 !             layertypes are explained below 
    96 !             icefrac = fraction of ice in icelayer
    97 !
    98 !     OUTPUTS: z = grid coordinates
    99 !              vectors ti and rhocv
    100 !***********************************************************************
    101   implicit none
    102   integer, intent(IN) :: nz, layertype
    103   real(8), intent(IN) :: zdepth, thIn, rhoc, porosity, icefrac
    104   real(8), intent(INOUT) :: z(nz)
    105   real(8), intent(OUT) :: ti(nz), rhocv(nz)
    106   integer j, b
    107   real(8) stretch, newrhoc, newti
    108   real(8), parameter :: NULL=0.
    109  
    110   if (zdepth>0 .and. zdepth<z(nz)) then
    111 
    112      select case (layertype)
    113      case (1)  ! interstitial ice
    114         call soilthprop(porosity,1.d0,rhoc,thIn,1,newrhoc,newti,NULL)
    115      case (2)  ! mostly ice (massive ice)
    116         call soilthprop(porosity,NULL,rhoc,thIn,2,newrhoc,newti,icefrac)
    117      case (3)  ! all ice (pure ice)
    118         call soilthprop(NULL,NULL,NULL,NULL,3,newrhoc,newti,NULL)
    119      case (4)  ! ice + rock + nothing else (ice-cemented soil)
    120         call soilthprop(porosity,NULL,rhoc,NULL,4,newrhoc,newti,NULL)
    121      case default
    122         error stop 'invalid layer type'
    123      end select
    124 
    125      ! Thermal skin depth is proportional to sqrt(kappa)
    126      ! thermal diffusivity kappa = k/(rho*c) = I^2/(rho*c)**2
    127      stretch = (newti/thIn)*(rhoc/newrhoc) ! ratio of sqrt(thermal diffusivity)
    128      
    129      b = 1
    130      do j=1,nz
    131         if (z(j)<=zdepth) then
    132            b = j+1
    133            rhocv(j) = rhoc
    134            ti(j) = thIn
    135         else
    136            rhocv(j) = newrhoc
    137            ti(j) = newti
    138         endif
    139         ! print *,j,z(j),ti(j),rhocv(j)
    140      enddo
    141      do j=b+1,nz
    142         z(j) = z(b) + stretch*(z(j)-z(b))
    143      enddo
    144      
    145      ! print *,'zmax=',z(nz),' b=',b,' stretch=',stretch
    146      ! print *,'depth at b-1, b ',z(b-1),z(b)
    147      ! print *,'ti1=',thIn,' ti2=',newti
    148      ! print *,'rhoc1=',rhoc,' rhoc2=',newrhoc
    149   endif
    150  
    151 end subroutine smartgrid
    152 

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

    r3525 r3526  
    492492== 19/11/2024 == JBC
    493493Computation of <soil thermal intertia> and <H2O mass subsurface/surface exchange> according to the presence of subsurface ice provided by the (Norbert's) dynamic method + few cleanings.
     494
     495== 19/11/2024 == JBC
     496Removing unused or redundant Norbert Schorghofer's subroutines (follow-up of r3493) + cleaning and some modifications of related subroutines.

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

    r2961 r3526  
    22!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    33!!!
    4 !!! Purpose: The module contains all the mathematical subroutine used in the PEM
     4!!! Purpose: The module contains all the mathematical SUBROUTINE used in the PEM
    55!!!           
    66!!! Author: Adapted from Schorgofer MSIM (N.S, Icarus 2010), impletented here by LL
    77!!!
    88!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    9   implicit none
     9
     10implicit none
     11
     12!=======================================================================
    1013  contains
     14!=======================================================================
    1115
    12 
    13 subroutine deriv1(z,nz,y,y0,ybot,dzY)
     16SUBROUTINE deriv1(z,nz,y,y0,ybot,dzY)
    1417!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1518!!!
    … …  
    1922!!!
    2023!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    21         implicit none
     24! first derivative of a function y(z) on irregular grid
     25! upper boundary conditions: y(0) = y0
     26! lower boundary condition.: yp = ybottom
    2227
    23 ! first derivative of a function y(z) on irregular grid
    24 ! upper boundary conditions: y(0)=y0
    25 ! lower boundary condition.: yp = ybottom
    26   integer, intent(IN) :: nz ! number of layer
    27   real, intent(IN) :: z(nz) ! depth layer
    28   real, intent(IN) :: y(nz) ! function which needs to be derived
    29   real, intent(IN) :: y0,ybot ! boundary conditions
    30   real, intent(OUT) :: dzY(nz) ! derivative of y w.r.t depth
     28implicit none
     29
     30integer, intent(in) :: nz      ! number of layer
     31real,    intent(in) :: z(nz)   ! depth layer
     32real,    intent(in) :: y(nz)   ! function which needs to be derived
     33real,    intent(in) :: y0,ybot ! boundary conditions
     34real, intent(out) :: dzY(nz) ! derivative of y w.r.t depth
    3135! local
    32   integer :: j
    33   real :: hm,hp,c1,c2,c3
     36integer :: j
     37real    :: hm, hp, c1, c2, c3
    3438
    35   hp = z(2)-z(1)
    36   c1 = z(1)/(hp*z(2))
    37   c2 = 1/z(1) - 1/(z(2)-z(1))
    38   c3 = -hp/(z(1)*z(2))
    39   dzY(1) = c1*y(2) + c2*y(1) + c3*y0
    40   do j=2,nz-1
    41      hp = z(j+1)-z(j)
    42      hm = z(j)-z(j-1)
    43      c1 = +hm/(hp*(z(j+1)-z(j-1)))
    44      c2 = 1/hm - 1/hp
    45      c3 = -hp/(hm*(z(j+1)-z(j-1)))
    46      dzY(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
    47   enddo
    48   dzY(nz) = (ybot - y(nz-1))/(2.*(z(nz)-z(nz-1)))
    49  return
    50 end subroutine deriv1
     39hp = z(2) - z(1)
     40c1 = z(1)/(hp*z(2))
     41c2 = 1/z(1) - 1/(z(2) - z(1))
     42c3 = -hp/(z(1)*z(2))
     43dzY(1) = c1*y(2) + c2*y(1) + c3*y0
     44do j = 2,nz - 1
     45    hp = z(j + 1) - z(j)
     46    hm = z(j) - z(j - 1)
     47    c1 = +hm/(hp*(z(j + 1) - z(j - 1)))
     48    c2 = 1/hm - 1/hp
     49    c3 = -hp/(hm*(z(j + 1) - z(j - 1)))
     50    dzY(j) = c1*y(j + 1) + c2*y(j) + c3*y(j - 1)
     51enddo
     52dzY(nz) = (ybot - y(nz - 1))/(2.*(z(nz) - z(nz - 1)))
    5153
     54END SUBROUTINE deriv1
    5255
     56!=======================================================================
    5357
    54 subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2)
     58SUBROUTINE deriv2_simple(z,nz,y,y0,yNp1,yp2)
    5559!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    5660!!!
    … …  
    6064!!!
    6165!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     66! second derivative y_zz on irregular grid
     67! boundary conditions: y(0) = y0, y(nz + 1) = yNp1
    6268
    63   ! second derivative y_zz on irregular grid
    64   ! boundary conditions: y(0)=y0, y(nz+1)=yNp1
    65   implicit none
    66   integer, intent(IN) :: nz
    67   real, intent(IN) :: z(nz),y(nz),y0,yNp1
    68   real, intent(OUT) :: yp2(nz)
    69   integer j
    70   real hm,hp,c1,c2,c3
     69implicit none
    7170
    72   c1 = +2./((z(2)-z(1))*z(2))
    73   c2 = -2./((z(2)-z(1))*z(1))
    74   c3 = +2./(z(1)*z(2))
    75   yp2(1) = c1*y(2) + c2*y(1) + c3*y0
     71integer, intent(in) :: nz
     72real,    intent(in) :: z(nz), y(nz), y0, yNp1
     73real, intent(out) :: yp2(nz)
     74integer :: j
     75real    :: hm, hp, c1, c2, c3
    7676
    77   do j=2,nz-1
    78      hp = z(j+1)-z(j)
    79      hm = z(j)-z(j-1)
    80      c1 = +2./(hp*(z(j+1)-z(j-1)))
    81      c2 = -2./(hp*hm)
    82      c3 = +2./(hm*(z(j+1)-z(j-1)))
    83      yp2(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
    84   enddo
    85   yp2(nz) = (yNp1 - 2*y(nz) + y(nz-1))/(z(nz)-z(nz-1))**2
    86   return
    87 end subroutine deriv2_simple
     77c1 = +2./((z(2) - z(1))*z(2))
     78c2 = -2./((z(2) - z(1))*z(1))
     79c3 = +2./(z(1)*z(2))
     80yp2(1) = c1*y(2) + c2*y(1) + c3*y0
     81do j = 2,nz - 1
     82    hp = z(j + 1) - z(j)
     83    hm = z(j) - z(j - 1)
     84    c1 = +2./(hp*(z(j + 1) - z(j - 1)))
     85    c2 = -2./(hp*hm)
     86    c3 = +2./(hm*(z(j + 1) - z(j-1)))
     87    yp2(j) = c1*y(j + 1) + c2*y(j) + c3*y(j - 1)
     88enddo
     89yp2(nz) = (yNp1 - 2*y(nz) + y(nz - 1))/(z(nz) - z(nz - 1))**2
    8890
     91END SUBROUTINE deriv2_simple
    8992
     93!=======================================================================
    9094
    91 subroutine  deriv1_onesided(j,z,nz,y,dy_zj)
     95SUBROUTINE  deriv1_onesided(j,z,nz,y,dy_zj)
    9296!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    9397!!!
    94 !!! Purpose:   First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
     98!!! Purpose: First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
    9599!!!           
    96100!!! Author: N.S (raw copy/paste from MSIM)
    … …  
    99103
    100104  implicit none
    101   integer, intent(IN) :: nz,j
    102   real, intent(IN) :: z(nz),y(nz)
    103   real, intent(out) :: dy_zj
    104   real h1,h2,c1,c2,c3
    105   if (j<1 .or. j>nz-2) then
    106      dy_zj = -1.
    107   else
    108      h1 = z(j+1)-z(j)
    109      h2 = z(j+2)-z(j+1)
    110      c1 = -(2*h1+h2)/(h1*(h1+h2))
    111      c2 =  (h1+h2)/(h1*h2)
    112      c3 = -h1/(h2*(h1+h2))
    113      dy_zj = c1*y(j) + c2*y(j+1) + c3*y(j+2)
    114   endif
    115   return
    116 end subroutine deriv1_onesided
    117105
     106integer, intent(in) :: nz, j
     107real,    intent(in) :: z(nz), y(nz)
     108real, intent(out) :: dy_zj
     109real :: h1, h2, c1, c2, c3
    118110
    119 subroutine  colint(y,z,nz,i1,i2,integral)
     111if (j < 1 .or. j > nz - 2) then
     112    dy_zj = -1.
     113else
     114    h1 = z(j + 1) - z(j)
     115    h2 = z(j + 2)- z(j + 1)
     116    c1 = -(2*h1 + h2)/(h1*(h1 + h2))
     117    c2 = (h1 + h2)/(h1*h2)
     118    c3 = -h1/(h2*(h1 + h2))
     119    dy_zj = c1*y(j) + c2*y(j + 1) + c3*y(j + 2)
     120endif
     121
     122END SUBROUTINE deriv1_onesided
     123
     124!=======================================================================
     125
     126PURE SUBROUTINE colint(y,z,nz,i1,i2,integral)
    120127!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    121128!!!
    … …  
    126133!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    127134 
    128   implicit none
    129   integer, intent(IN) :: nz, i1, i2
    130   real, intent(IN) :: y(nz), z(nz)
    131   real,intent(out) :: integral
    132   integer i
    133   real dz(nz)
     135implicit none
    134136 
    135   dz(1) = (z(2)-0.)/2
    136   do i=2,nz-1
    137      dz(i) = (z(i+1)-z(i-1))/2.
    138   enddo
    139   dz(nz) = z(nz)-z(nz-1)
    140   integral = sum(y(i1:i2)*dz(i1:i2))
    141 end subroutine colint
     137integer, intent(in) :: nz, i1, i2
     138real,    intent(in) :: y(nz), z(nz)
     139real, intent(out) :: integral
     140integer i
     141real dz(nz)
    142142
     143dz(1) = (z(2) - 0.)/2
     144do i = 2,nz - 1
     145    dz(i) = (z(i + 1) - z(i - 1))/2.
     146enddo
     147dz(nz) = z(nz) - z(nz - 1)
     148integral = sum(y(i1:i2)*dz(i1:i2))
    143149
     150END SUBROUTINE colint
    144151
    145 
     152!=======================================================================
    146153
    147154SUBROUTINE findroot(y1,y2,z1,z2,zr)
    148         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    149         !!!
    150         !!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
    151         !!!
    152         !!! Author: LL
    153         !!!
    154         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     155!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     156!!!
     157!!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
     158!!!
     159!!! Author: LL
     160!!!
     161!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    155162
    156         implicit none
    157         real,intent(in) :: y1,y2 ! difference between surface water density and at depth  [kg/m^3]
    158         real,intent(in) :: z1,z2 ! depth [m]
    159         real,intent(out) :: zr ! depth at which we have zero
    160         zr = (y1*z2 - y2*z1)/(y1-y2)
    161         RETURN
    162         end subroutine findroot
     163implicit none
     164
     165real, intent(in) :: y1, y2 ! difference between surface water density and at depth  [kg/m^3]
     166real, intent(in) :: z1, z2 ! depth [m]
     167real, intent(out) :: zr    ! depth at which we have zero
     168
     169zr = (y1*z2 - y2*z1)/(y1 - y2)
     170
     171END SUBROUTINE findroot
    163172
    164173end module math_mod
Note: See TracChangeset for help on using the changeset viewer.