Changeset 3526 for trunk/LMDZ.COMMON/libf/evolution
- Timestamp:
- Nov 19, 2024, 5:54:18 PM (5 weeks ago)
- 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 2 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3 2 !!! … … 171 170 172 171 end subroutine dyn_ss_ice_m 173 174 175 -
trunk/LMDZ.COMMON/libf/evolution/NS_fast_modules.F90
r3493 r3526 1 2 1 module allinterfaces 3 2 ! interfaces from Fortran 90 subroutines and functions … … 86 85 real(8), intent(OUT) :: ypp(nz), zdepthG 87 86 integer, intent(INOUT) :: typeG 88 real(8), external :: zint , deriv1_onesided, colint87 real(8), external :: zint 89 88 end subroutine depths_avmeth 90 89 end interface … … 96 95 real(8) constriction 97 96 end function constriction 98 end interface99 100 interface101 subroutine assignthermalproperties(nz,thIn,rhoc,ti,rhocv, &102 & typeT,icefrac,porosity,porefill)103 implicit none104 integer, intent(IN) :: nz105 integer, intent(IN), optional :: typeT106 real(8), intent(IN), optional :: icefrac107 real(8), intent(IN) :: thIn, rhoc108 real(8), intent(IN), optional :: porosity, porefill(nz)109 real(8), intent(OUT) :: ti(nz), rhocv(nz)110 end subroutine assignthermalproperties111 end interface112 113 interface114 pure subroutine icechanges_poreonly(nz,z,typeF,typeG,avdrhoP,ypp,B,porefill)115 implicit none116 integer, intent(IN) :: nz, typeF, typeG117 real(8), intent(IN) :: z(nz), ypp(nz), avdrhoP, B118 real(8), intent(INOUT) :: porefill(nz)119 end subroutine icechanges_poreonly120 97 end interface 121 98 … … 131 108 end interface 132 109 133 interface134 subroutine compactoutput(unit,porefill,nz)135 implicit none136 integer, intent(IN) :: unit,nz137 real(8), intent(IN) :: porefill(nz)138 end subroutine compactoutput139 end interface140 141 110 !end of fast_subs_univ 142 !begin derivs.f90143 144 interface145 subroutine deriv1(z,nz,y,y0,yNp1,yp)146 implicit none147 integer, intent(IN) :: nz148 real(8), intent(IN) :: z(nz),y(nz),y0,yNp1149 real(8), intent(OUT) :: yp(nz)150 end subroutine deriv1151 end interface152 153 interface154 subroutine deriv2_full(z,nz,a,b,a0,b0,bNp1,yp2)155 implicit none156 integer, intent(IN) :: nz157 real(8), intent(IN) :: z(nz),a(nz),b(nz),a0,b0,bNp1158 real(8), intent(OUT) :: yp2(nz)159 end subroutine deriv2_full160 end interface161 162 interface163 subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2)164 implicit none165 integer, intent(IN) :: nz166 real(8), intent(IN) :: z(nz),y(nz),y0,yNp1167 real(8), intent(OUT) :: yp2(nz)168 end subroutine deriv2_simple169 end interface170 interface171 real(8) pure function deriv1_onesided(j,z,nz,y)172 implicit none173 integer, intent(IN) :: nz,j174 real(8), intent(IN) :: z(nz),y(nz)175 end function deriv1_onesided176 end interface177 178 !end of derivs.f90179 !begin fast_subs_exper.f90180 181 interface182 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, icedensity185 implicit none186 integer, intent(IN) :: nz187 real(8), intent(IN) :: bigstep188 real(8), intent(IN) :: thIn, rhoc, z(NMAX), porosity, pfrost189 real(8), intent(INOUT) :: zdepthT, zdepthF, porefill(nz)190 real(8), intent(OUT) :: zdepthE191 real(8), intent(IN) :: icefrac, Diff, Tmean, Tampl192 real(8), intent(OUT) :: zdepthG193 end subroutine icelayer_exper194 end interface195 196 interface197 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, marsDay201 implicit none202 integer, intent(IN) :: nz, typeT203 real(8), intent(IN) :: z(NMAX), ti(NMAX), rhocv(NMAX), pfrost204 real(8), intent(IN) :: Tmean, Tampl205 real(8), intent(OUT) :: avdrho, avdrhoP206 real(8), intent(OUT) :: zdepthE207 integer, intent(OUT) :: typeF208 real(8), intent(INOUT) :: zdepthF209 real(8), intent(OUT) :: ypp(nz)210 real(8), intent(IN) :: porefill(nz)211 real(8), intent(IN) :: B212 integer, intent(OUT) :: typeG213 real(8), intent(OUT) :: zdepthG214 real(8), external :: Tsurface, psv215 real(8), external :: equildepth216 end subroutine ajsub_exper217 end interface218 219 !end fast_subs_exper220 221 111 ! Other 222 interface223 pure function flux_mars77(R,decl,HA,albedo,fracir,fracscat)224 implicit none225 real*8 flux_mars77226 real*8, intent(IN) :: R,decl,HA,albedo,fracIR,fracScat227 end function flux_mars77228 end interface229 230 interface231 pure function colint(y,z,nz,i1,i2)232 implicit none233 integer, intent(IN) :: nz, i1, i2234 real(8), intent(IN) :: y(nz),z(nz)235 real(8) colint236 end function colint237 end interface238 239 112 interface 240 113 subroutine dyn_ss_ice_m(ssi_depth_in,T1,T_in,nsoil,thIn,p0,pfrost,porefill_in,porefill,ssi_depth) … … 250 123 end interface 251 124 252 interface253 subroutine dyn_ss_ice_m_wrapper(ngrid,nsoil,tHIn,p0,pfrost,T_in,ssi_depth_in,porefill_in,porefill,ssi_depth)254 implicit none255 integer, intent(IN) :: nsoil,ngrid256 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_wrapper263 end interface264 265 125 end module allinterfaces -
trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90
r3512 r3526 294 294 295 295 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 296 real*8 function tfrostco2(p) 297 ! the inverse of function psvco2 298 ! input is pressure in Pascal, output is temperature in Kelvin 299 implicit none 300 real*8 p 301 tfrostco2 = 3182.48/(23.3494+log(100./p)) 302 end function -
trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90
r3493 r3526 55 55 !*********************************************************************** 56 56 use allinterfaces, except_this_one => depths_avmeth 57 use math_mod, only: deriv2_simple, deriv1_onesided, deriv1, colint 57 58 implicit none 58 59 integer, intent(IN) :: nz, typeT … … 118 119 call deriv1(z,nz,eta(:),1.d0,eta(nz-1),ap) 119 120 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) 121 122 ! print *,typeP,ap(typeP),ap_one 122 123 ap(typeP)=ap_one … … 151 152 endif 152 153 if (typeG>0 .and. typeT<0) then 153 c umfillabove = colint(porefill(:)/eta(:),z,nz,typeG-1,nz)154 call colint(porefill(:)/eta(:),z,nz,typeG-1,nz,cumfillabove) 154 155 newtypeG = -9 155 156 do i=typeG,nz 156 157 if (minval(eta(i:nz))<=0.) cycle ! eta=0 means completely full 157 c umfill=colint(porefill(:)/eta(:),z,nz,i,nz)158 call colint(porefill(:)/eta(:),z,nz,i,nz,cumfill) 158 159 if (cumfill<yp(i)*18./8314.*B) then ! usually executes on i=typeG 159 160 if (i>typeG) then … … 194 195 195 196 196 pure subroutine icechanges_poreonly(nz,z,typeF,typeG,avdrhoP,ypp,B,porefill)197 use allinterfaces, except_this_one => icechanges_poreonly198 implicit none199 integer, intent(IN) :: nz, typeF, typeG200 real(8), intent(IN) :: z(nz), ypp(nz), avdrhoP, B201 real(8), intent(INOUT) :: porefill(nz)202 integer j, erase, newtypeP, ub203 real(8) integ204 205 !----retreat206 ! avdrhoP>0 is outward loss from zdepthP207 ! avdrhoP<0 means gain at zdepthP or no ice anywhere208 if (avdrhoP>0.) then209 erase=0210 do j=1,nz211 if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF212 integ = colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j)213 erase = j214 if (integ>B*avdrhoP*18./8314.) exit215 end do216 if (erase>0) porefill(1:erase)=0.217 endif218 219 ! new depth220 newtypeP = -9221 do j=1,nz222 if (porefill(j)>0.) then223 newtypeP = j ! first point with ice224 exit225 endif226 enddo227 228 !----diffusive filling229 ub = typeF230 if (newtypeP>0 .and. typeF>0 .and. newtypeP<ub) ub=newtypeP231 if (ub>0) then232 do j=ub,nz233 ! B=Diff/(porosity*icedensity)*86400*365.24*bigstep234 porefill(j) = porefill(j) + B*ypp(j)235 if (porefill(j)<0.) porefill(j)=0.236 if (porefill(j)>1.) porefill(j)=1.237 enddo238 end if239 240 !----enact bottom boundary241 if (typeG>0) porefill(typeG:nz)=0.242 243 end subroutine icechanges_poreonly244 245 246 247 197 pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, & 248 198 & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG) … … 254 204 use miscparameters, only : icedensity 255 205 use allinterfaces, except_this_one => icechanges 206 use math_mod, only: colint 256 207 implicit none 257 208 integer, intent(IN) :: nz, typeF, typeG … … 293 244 if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF 294 245 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) 296 247 erase = j 297 248 if (integ>B*avdrhoP*18./8314.) exit … … 332 283 end if 333 284 end subroutine icechanges 334 335 336 subroutine assignthermalproperties(nz,thIn,rhoc, &337 & ti,rhocv,typeT,icefrac,porosity,porefill)338 !*********************************************************339 ! assign thermal properties of soil340 !*********************************************************341 implicit none342 integer, intent(IN) :: nz343 integer, intent(IN), optional :: typeT344 real(8), intent(IN), optional :: icefrac345 real(8), intent(IN) :: thIn, rhoc346 real(8), intent(IN), optional :: porosity, porefill(nz)347 real(8), intent(OUT) :: ti(nz), rhocv(nz)348 integer j349 real(8) newrhoc, newti, fill350 real(8), parameter :: NULL=0.351 352 ti(1:nz) = thIn353 rhocv(1:nz) = rhoc354 if (typeT>0) then355 call soilthprop(porosity,NULL,rhoc,thIn,2,newrhoc,newti,icefrac)356 rhocv(typeT:nz) = newrhoc357 ti(typeT:nz) = newti358 endif359 do j=1,nz360 fill = porefill(j) ! off by half point361 if (fill>0. .and. (typeT<0 .or. (typeT>0 .and. j<typeT))) then362 call soilthprop(porosity,fill,rhoc,thIn,1,rhocv(j),ti(j),NULL)363 endif364 enddo365 end subroutine assignthermalproperties366 367 368 369 subroutine compactoutput(unit,porefill,nz)370 implicit none371 integer, intent(IN) :: unit,nz372 real(8), intent(IN) :: porefill(nz)373 integer j374 do j=1,nz375 if (porefill(j)==0.) then376 write(unit,'(1x,f2.0)',advance='no') porefill(j)377 else378 write(unit,'(1x,f7.5)',advance='no') porefill(j)379 endif380 enddo381 write(unit,"('')")382 end subroutine compactoutput383 -
trunk/LMDZ.COMMON/libf/evolution/NS_misc_params.F90
r3493 r3526 12 12 ! length of Earth year in days = 365.24 13 13 end module miscparameters 14 -
trunk/LMDZ.COMMON/libf/evolution/NS_soilthprop_mars.F90
r3493 r3526 76 76 77 77 end 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 appropriate84 ! values of thermal inertia ti and rho*c in icy layer85 !86 ! INPUTS:87 ! nz = number of grid points88 ! z = grid spacing for dry regolith89 ! (will be partially overwritten)90 ! zdepth = depth where ice table starts91 ! negative values indicate no ice92 ! 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 volume95 ! layertypes are explained below96 ! icefrac = fraction of ice in icelayer97 !98 ! OUTPUTS: z = grid coordinates99 ! vectors ti and rhocv100 !***********************************************************************101 implicit none102 integer, intent(IN) :: nz, layertype103 real(8), intent(IN) :: zdepth, thIn, rhoc, porosity, icefrac104 real(8), intent(INOUT) :: z(nz)105 real(8), intent(OUT) :: ti(nz), rhocv(nz)106 integer j, b107 real(8) stretch, newrhoc, newti108 real(8), parameter :: NULL=0.109 110 if (zdepth>0 .and. zdepth<z(nz)) then111 112 select case (layertype)113 case (1) ! interstitial ice114 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 default122 error stop 'invalid layer type'123 end select124 125 ! Thermal skin depth is proportional to sqrt(kappa)126 ! thermal diffusivity kappa = k/(rho*c) = I^2/(rho*c)**2127 stretch = (newti/thIn)*(rhoc/newrhoc) ! ratio of sqrt(thermal diffusivity)128 129 b = 1130 do j=1,nz131 if (z(j)<=zdepth) then132 b = j+1133 rhocv(j) = rhoc134 ti(j) = thIn135 else136 rhocv(j) = newrhoc137 ti(j) = newti138 endif139 ! print *,j,z(j),ti(j),rhocv(j)140 enddo141 do j=b+1,nz142 z(j) = z(b) + stretch*(z(j)-z(b))143 enddo144 145 ! print *,'zmax=',z(nz),' b=',b,' stretch=',stretch146 ! print *,'depth at b-1, b ',z(b-1),z(b)147 ! print *,'ti1=',thIn,' ti2=',newti148 ! print *,'rhoc1=',rhoc,' rhoc2=',newrhoc149 endif150 151 end subroutine smartgrid152 -
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3525 r3526 492 492 == 19/11/2024 == JBC 493 493 Computation 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 496 Removing 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 2 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3 3 !!! 4 !!! Purpose: The module contains all the mathematical subroutineused in the PEM4 !!! Purpose: The module contains all the mathematical SUBROUTINE used in the PEM 5 5 !!! 6 6 !!! Author: Adapted from Schorgofer MSIM (N.S, Icarus 2010), impletented here by LL 7 7 !!! 8 8 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 9 implicit none 9 10 implicit none 11 12 !======================================================================= 10 13 contains 14 !======================================================================= 11 15 12 13 subroutine deriv1(z,nz,y,y0,ybot,dzY) 16 SUBROUTINE deriv1(z,nz,y,y0,ybot,dzY) 14 17 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 15 18 !!! … … 19 22 !!! 20 23 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 22 27 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 28 implicit none 29 30 integer, intent(in) :: nz ! number of layer 31 real, intent(in) :: z(nz) ! depth layer 32 real, intent(in) :: y(nz) ! function which needs to be derived 33 real, intent(in) :: y0,ybot ! boundary conditions 34 real, intent(out) :: dzY(nz) ! derivative of y w.r.t depth 31 35 ! local 32 33 real :: hm,hp,c1,c2,c336 integer :: j 37 real :: hm, hp, c1, c2, c3 34 38 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 39 hp = z(2) - z(1) 40 c1 = z(1)/(hp*z(2)) 41 c2 = 1/z(1) - 1/(z(2) - z(1)) 42 c3 = -hp/(z(1)*z(2)) 43 dzY(1) = c1*y(2) + c2*y(1) + c3*y0 44 do 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) 51 enddo 52 dzY(nz) = (ybot - y(nz - 1))/(2.*(z(nz) - z(nz - 1))) 51 53 54 END SUBROUTINE deriv1 52 55 56 !======================================================================= 53 57 54 subroutinederiv2_simple(z,nz,y,y0,yNp1,yp2)58 SUBROUTINE deriv2_simple(z,nz,y,y0,yNp1,yp2) 55 59 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 56 60 !!! … … 60 64 !!! 61 65 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 66 ! second derivative y_zz on irregular grid 67 ! boundary conditions: y(0) = y0, y(nz + 1) = yNp1 62 68 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 69 implicit none 71 70 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 71 integer, intent(in) :: nz 72 real, intent(in) :: z(nz), y(nz), y0, yNp1 73 real, intent(out) :: yp2(nz) 74 integer :: j 75 real :: hm, hp, c1, c2, c3 76 76 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 77 c1 = +2./((z(2) - z(1))*z(2)) 78 c2 = -2./((z(2) - z(1))*z(1)) 79 c3 = +2./(z(1)*z(2)) 80 yp2(1) = c1*y(2) + c2*y(1) + c3*y0 81 do 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) 88 enddo 89 yp2(nz) = (yNp1 - 2*y(nz) + y(nz - 1))/(z(nz) - z(nz - 1))**2 88 90 91 END SUBROUTINE deriv2_simple 89 92 93 !======================================================================= 90 94 91 subroutinederiv1_onesided(j,z,nz,y,dy_zj)95 SUBROUTINE deriv1_onesided(j,z,nz,y,dy_zj) 92 96 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 93 97 !!! 94 !!! Purpose: 98 !!! Purpose: First derivative of function y(z) at z(j) one-sided derivative on irregular grid 95 99 !!! 96 100 !!! Author: N.S (raw copy/paste from MSIM) … … 99 103 100 104 implicit none 101 integer, intent(IN) :: nz,j102 real, intent(IN) :: z(nz),y(nz)103 real, intent(out) :: dy_zj104 real h1,h2,c1,c2,c3105 if (j<1 .or. j>nz-2) then106 dy_zj = -1.107 else108 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 endif115 return116 end subroutine deriv1_onesided117 105 106 integer, intent(in) :: nz, j 107 real, intent(in) :: z(nz), y(nz) 108 real, intent(out) :: dy_zj 109 real :: h1, h2, c1, c2, c3 118 110 119 subroutine colint(y,z,nz,i1,i2,integral) 111 if (j < 1 .or. j > nz - 2) then 112 dy_zj = -1. 113 else 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) 120 endif 121 122 END SUBROUTINE deriv1_onesided 123 124 !======================================================================= 125 126 PURE SUBROUTINE colint(y,z,nz,i1,i2,integral) 120 127 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 121 128 !!! … … 126 133 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 127 134 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) 135 implicit none 134 136 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 137 integer, intent(in) :: nz, i1, i2 138 real, intent(in) :: y(nz), z(nz) 139 real, intent(out) :: integral 140 integer i 141 real dz(nz) 142 142 143 dz(1) = (z(2) - 0.)/2 144 do i = 2,nz - 1 145 dz(i) = (z(i + 1) - z(i - 1))/2. 146 enddo 147 dz(nz) = z(nz) - z(nz - 1) 148 integral = sum(y(i1:i2)*dz(i1:i2)) 143 149 150 END SUBROUTINE colint 144 151 145 152 !======================================================================= 146 153 147 154 SUBROUTINE findroot(y1,y2,z1,z2,zr) 148 149 150 151 152 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 155 162 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 163 implicit none 164 165 real, intent(in) :: y1, y2 ! difference between surface water density and at depth [kg/m^3] 166 real, intent(in) :: z1, z2 ! depth [m] 167 real, intent(out) :: zr ! depth at which we have zero 168 169 zr = (y1*z2 - y2*z1)/(y1 - y2) 170 171 END SUBROUTINE findroot 163 172 164 173 end module math_mod
Note: See TracChangeset
for help on using the changeset viewer.