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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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.