#if defined O_mtlm_nsoil subroutine MTLM_STATE (DT, G, ETRAN_L, MELTHEAT) # if defined O_mtlm !----------------------------------------------------------------------- ! Routine to update land surface prognostic variables ! (soil moisture, soil temperature, lying snow mass). ! Also diagnoses runoff and snowmelt !********************************************************************** ! this file is based on code that may have had the following copyright: ! (c) CROWN COPYRIGHT 1997, U.K. METEOROLOGICAL OFFICE. ! Permission has been granted by the authors to the public to copy ! and use this software without charge, provided that this Notice and ! any statement of authorship are reproduced on all copies. Neither the ! Crown nor the U.K. Meteorological Office makes any warranty, express ! or implied, or assumes any liability or responsibility for the use of ! this software. !********************************************************************** !----------------------------------------------------------------------- implicit none include "size.h" include "mtlm.h" include "mtlm_data.h" !~ Loop/spatial variables ~ ! I, L, N, X = Loop counters ! Z = Depth counter ! NPERMA = Soil layer containing permafrost integer I, L, N, NPERMA, X real Z !~ Constants ~ ! DT = Timestep (s). ! FILL = Fill value for permafrost array ! SECSDAY = Seconds per day ! SLOPEARRAY= Value of L_slope.nc to use for wetlands diagnosis ! WETVAL = Value of SU that wet areas must exceed to be wetlands real DT, WETVAL real FILL, SECSDAY integer SLOPEARRAY PARAMETER (FILL = -9.9692099683868690e+36, SLOPEARRAY = 4) PARAMETER (SECSDAY = 86400, WETVAL = 0.65) !~ Other variables ~ ! DZ_EFF = Effective layer thickness (excess ice) (m) ! ETRAN_L = Evapotranspiration from each layer (kg/m2/s). ! G = Ground heat flux (W/m2). ! HCAP_A = Actual heat capacity (dry soil + liquid moisture) ! DRAINAGE = Flow of water through ground water (kg/m2/s) ! INFIL = Maximum rate of infiltration into soil ! MELTABLE = Amount of lying_snow it is possbile to melt (kg/m2) ! MELTHEAT = Heat available to melt snow (W/m2) ! W_FLUX = Water flux (kg/m2/s) into layer N real ETRAN_L(POINTS, NGND), G(POINTS), INFIL(POINTS) real DRAINAGE(POINTS) real W_FLUX(POINTS, 0:NGND), HCAP_A real MELTABLE, MELTHEAT(POINTS), DZ_EFF(POINTS, NGND) real SUBSOIL_L(POINTS, NGND), ADDMELT real DELTA, MAXFROZEN(POINTS) DZ_EFF(:,:) = 0 ! These 4 lines of code went missing when nsoil was integrated into model by ! Eby in autumn 2013. Need to be here for diagnostics and Kovin diffusion to ! work properly PERMA(:) = FILL PFTHICK(:) = 0. ACTIVELT(:) = 0. TALIK(:) = 0. do I=1,LAND_PTS L = LAND_INDEX(I) # if defined O_mtlm_exstruct do N = 1, NSOIL(L) DZ_EFF(L,N) = DZ_GND(N) + LAYERICE(L,N)/RHO_I enddo do N = NSOIL(L)+1, NGND DZ_EFF(L,N) = DZ_GND(N) enddo # else DZ_EFF(L,:) = DZ_GND(:) # endif !------------------------------------------------------------------------- ! Diagnose snow melt and update the lying snow ! Order of priority for snow fluxes: ! i) Sublimation ! ii) Surface melt ! iii) Subsurface melt !------------------------------------------------------------------------- ! i) Sublimation + snowfall fluxes SNOWMELT(L) = 0.0 LYING_SNOW(L) = LYING_SNOW(L) + (SNOW(L) - SUBSNOW(L))*DT if (LYING_SNOW(L) .lt. 0.) then ! If negative snow, change extra sublimation to evaporation SUBSNOW(L) = SUBSNOW(L) + LYING_SNOW(L)/DT ETRAN_L(L,1) = ETRAN_L(L,1) - LYING_SNOW(L)/DT G(L) = G(L) + LHF*LYING_SNOW(L) / DT LYING_SNOW(L) = 0. MELTABLE = 0. else ! Still snow left to melt ! ii) Surface melt MELTABLE = 0. MELTABLE = min(MELTHEAT(L)*DT/LHF, LYING_SNOW(L)) LYING_SNOW(L) = LYING_SNOW(L) - MELTABLE SNOWMELT(L) = SNOWMELT(L) + MELTABLE / DT endif G(L) = G(L) + MELTHEAT(L)- (MELTABLE * LHF / DT) ! Adjust G in case too much melt heat ! iii) Subsurface melt if ((LYING_SNOW(L).gt.0.0).and.(TGND(L,1).gt.ZERODEGC)) then HCAP_A = HCAP_D(L,1) + HCAP_W * RHO_W * V_SAT(L,1) * SU(L,1) if ((HCAP_A*DZ_GND(1)*(TGND(L,1)-ZERODEGC)/LHF) & .gt. LYING_SNOW(L)) then TGND(L,1) = TGND(L,1) - (LYING_SNOW(L)*LHF) & / (HCAP_A*DZ_GND(1)) SNOWMELT(L) = SNOWMELT(L) + LYING_SNOW(L)/DT LYING_SNOW(L) = 0.! else SNOWMELT(L) = SNOWMELT(L) + HCAP_A*DZ_GND(1) & * (TGND(L,1)-ZERODEGC)/(LHF*DT) LYING_SNOW(L) = LYING_SNOW(L) - HCAP_A*DZ_GND(1) & * (TGND(L,1) - ZERODEGC) / LHF TGND(L,1) = ZERODEGC endif endif SNOW_HT(L) = LYING_SNOW(L)/RHO_S SURF_ROFF(L) = 0. W_FLUX(L,0) = RAIN(L) + SNOWMELT(L) enddo !------------------------------------------------------------------- ! Call the soil hydrology code to diagnose water fluxes between soil ! layers and through the bottom of the soil column !------------------------------------------------------------------- call SOIL_HYD (ETRAN_L, DT, DRAINAGE, W_FLUX, SUBSOIL_L) !------------------------------------------------------------------- ! Call the soil temperature code to diagnose changes in soil layer ! temperatures, heat fluxes between layers and phase changes of water !------------------------------------------------------------------- call SOIL_TEMP (DT, G, DRAINAGE, W_FLUX, ETRAN_L, SUBSOIL_L) do I = 1, LAND_PTS L = LAND_INDEX(I) ! Add drainage to surface runoff since the ocean can't tell the difference SURF_ROFF(L) = SURF_ROFF(L) + DRAINAGE(L) !----------------------------------------------------------------- ! Update permafrost diagnostics !----------------------------------------------------------------- ! Update the number of days that the ground has been frozen / unfrozen do N=1,NGND if (TGND(L,N) .lt. 273.15) then DFROZEN_C(L,N) = DFROZEN_C(L,N) + DT/SECSDAY DFROZEN_A(L,N) = DFROZEN_A(L,N) + DT/SECSDAY else DFROZEN_C(L,N) = 0. endif if (TGND(L,N) .gt. 273.15) then DUNFROZEN_C(L,N) = DUNFROZEN_C(L,N) + DT/SECSDAY else DUNFROZEN_C(L,N) = 0. endif enddo Z = 0. NPERMA = 0. do N=1,NGND if (DFROZEN_C(L,N) .gt. 730.) then ! Have hit permafrost in the current layer PERMA(L) = Z NPERMA = N Z = Z + DZ_EFF(L,N) exit else ! Otherwise, advance to next layer Z = Z + DZ_EFF(L,N) endif enddo if (NPERMA .lt. NGND) then ! Have found permafrost, but above bottommost layer do N=NPERMA+1, NGND if (DFROZEN_C(L,N) .lt. 730.) then PFTHICK(L) = Z-PERMA(L) exit else Z = Z + DZ_EFF(L,N) endif if (N .eq. NGND) then PFTHICK(L) = Z - PERMA(L) exit endif enddo elseif (DFROZEN_C(L,NGND) .gt. 730) then ! Permafrost just in the bottom layer PFTHICK(L) = DZ_EFF(L,NGND) endif ! find active layer thickness and talik thickness (for taliks above first PF layer) if (NPERMA .gt. 1) then ! Must have found permafrost, but there will be an ALT do N=1,NPERMA-1 if (DFROZEN_C(L,N) .lt. 730. .and. & DUNFROZEN_C(L,N) .lt. 730.) then ACTIVELT(L) = ACTIVELT(L) + DZ_EFF(L,N) else exit endif enddo do N=1,NPERMA-1 if (DUNFROZEN_C(L,N) .gt. 730.) then TALIK(L) = TALIK(L) + DZ_EFF(L,N) endif enddo endif # if defined O_mtlm_exslab if (MSLAB(L) .gt. 0.) then if (NPERMA .lt. NSOIL(L)+2) then ! PF found starting above the slab or just below slab PFTHICK(L) = PFTHICK(L) + MSLAB(L)/RHO_I endif endif # endif enddo !------------------------------------------------------------------ ! Diagnose wetlands !------------------------------------------------------------------ YEARCOUNT = YEARCOUNT + DT/SECSDAY do I = 1, LAND_PTS L = LAND_INDEX(I) if (SU(L,1) .ge. WETVAL) then DWET(L) = DWET(L) + DT/SECSDAY endif if (LYING_SNOW(L) .gt. 0.) then DSNOW(L) = DSNOW(L) + DT/SECSDAY endif enddo if ((YEARCOUNT .ge. 365.)) then WETMAP(:,:) = 0. MAXFROZEN = MAXVAL(DFROZEN_A,2) FROZENMAP(:,:) = 0. do I = 1,LAND_PTS L = LAND_INDEX(I) DWET(L) = DWET(L) * 365./YEARCOUNT if ((DWET(L) .gt. 0).and.(DWET(L).le. 365)) then ANN_DWET(L) = DWET(L) else ANN_DWET(L) = 0. endif if (DWET(L) .gt. 0.) then WETMAP(L,1) = SLOPEMAP(L,SLOPEARRAY) if (DWET(L) .lt. 360.) then WETMAP(L,3) = SLOPEMAP(L,SLOPEARRAY) else WETMAP(L,2) = SLOPEMAP(L,SLOPEARRAY) endif endif if (PERMA(L) .gt. 0.) then FROZENMAP(L,1) = 1. elseif (MAXFROZEN(L) .gt. 15.) then FROZENMAP(L,2) = 1. elseif (MAXFROZEN(L) .ge. 1.) then FROZENMAP(L,3) = 1. elseif (MAXFROZEN(L) .gt. 0.) then FROZENMAP(L,4) = 1. endif if (DSNOW(L) .gt. 15.) then FROZENMAP(L,5) = 1. elseif (DSNOW(L) .ge. 1.) then FROZENMAP(L,6) = 1. elseif (DSNOW(L). gt. 0.) then FROZENMAP(L,7) = 1. endif DWET(L) = 0. DFROZEN_A(L,:) = 0. DSNOW(L) = 0. enddo # if defined O_mtlm_kovdiff YEARKOV_SW = 1 # endif YEARCOUNT = 0. endif !------------------------------------------------------------------ ! Update Schaefer Dmin !------------------------------------------------------------------- do I = 1, LAND_PTS L = LAND_INDEX(I) if(INT(SCHDMIN(L)) .lt. SCHDMAX) then do N = 1,SCHDMAX if(INT(SCHDMIN(L)) .lt. N) then if(TGND(L,N) .gt. 273.15) then SCHDMIN(L)=REAL(N) endif endif enddo endif enddo !------------------------------------------------------------------- ! Reboot Schaefer Dmin during spinup b/c memory quantity !------------------------------------------------------------------- # if defined O_mtlm_schaef_spin if (ISTEP .eq. STEP_DAY) then L = SCH_COUNT/STEP_DAY if (MOD(L,SCH_REST) .eq. 0.) then SCHDMIN(:)=0. SCH_COUNT = 0 write(*,*) 'SCHDMIN has been reboted' endif endif # endif # endif return end #else subroutine MTLM_STATE (POINTS, LAND_PTS, LAND_INDEX, DZ_SOIL &, HCAP_SOIL, KS, THETA_SAT, LF, TM, TIMESTEP &, G, RAIN, SNOW, E, ESUB, M, LYING_SNOW, TS1 &, RUNOFF, SNOWMELT, MNEG) # if defined O_mtlm !----------------------------------------------------------------------- ! Routine to update land surface prognostic variables ! (soil moisture, soil temperature, lying snow mass). ! Also diagnoses runoff and snowmelt !********************************************************************** ! this file is based on code that may have had the following copyright: ! (c) CROWN COPYRIGHT 1997, U.K. METEOROLOGICAL OFFICE. ! Permission has been granted by the authors to the public to copy ! and use this software without charge, provided that this Notice and ! any statement of authorship are reproduced on all copies. Neither the ! Crown nor the U.K. Meteorological Office makes any warranty, express ! or implied, or assumes any liability or responsibility for the use of ! this software. !********************************************************************** !----------------------------------------------------------------------- implicit none ! POINTS = IN Maximum number of land points. ! LAND_PTS = IN Number of land points. ! LAND_INDEX = IN Indices of landpoints integer POINTS, LAND_PTS, LAND_INDEX(POINTS) integer I, L ! Surface parameters ! DZ_SOIL = IN Soil layer thickness (m). ! HCAP_SOIL = IN Soil heat capacity (W/m3/K). ! TIMESTEP = IN Timestep (s). ! KS = IN Saturated hydraulic conductivity ! THETA_SAT = IN Saturated volumetric soil moisture concentration ! (m3/m3). ! LF = Latent heat of fusion (J/kg). ! TM = Melting point of fresh water (K). ! G = IN Ground heat flux (W/m2). real DZ_SOIL, HCAP_SOIL, TIMESTEP, KS, THETA_SAT, LF, TM real G(POINTS) ! Driving variables ! RAIN = IN Rainfall (kg/m2/s). ! SNOW = IN Snowfall (kg/m2/s). ! E = IN Evapotranspiration (kg/m2/s). ! ESUB = IN Sublimation (kg/m2/s). ! M = INOUT Soil moisture (kg/m2). ! NEG = INOUT Negative Soil moisture (kg/m2). ! LYING_SNOW = INOUT Lying snow (kg/m2). ! TS1 = INOUT Soil temperature (K). ! RUNOFF = OUT Runoff (kg/m2/s). ! SNOWMELT = OUT Snow melt (kg/m2/s). real RAIN(POINTS), SNOW(POINTS), E(POINTS), ESUB(POINTS) real M(POINTS), MNEG(POINTS), LYING_SNOW(POINTS) real TS1(POINTS), RUNOFF(POINTS), SNOWMELT(POINTS) ! Local parameters ! B = Clapp-Hornberger exponent. real B parameter (B=6.6) !---------------------------------------------------------------------- ! Update the soil temperature and diagnose snowmelt. !---------------------------------------------------------------------- !CDIR NODEP do I=1,LAND_PTS L = LAND_INDEX(I) TS1(L) = TS1(L) + TIMESTEP* G(L) / (DZ_SOIL*HCAP_SOIL) SNOWMELT(L) = 0.0 if ((LYING_SNOW(L).gt.0.0) .and. (TS1(L).gt.TM)) then SNOWMELT(L) = HCAP_SOIL*DZ_SOIL*(TS1(L)-TM) & / (LF*TIMESTEP) if ((SNOWMELT(L)).gt.(SNOW(L)-ESUB(L)+LYING_SNOW(L)/TIMESTEP)) & then ! limit snowmelt to amount of snow and fix soil temperature SNOWMELT(L) = SNOW(L)-ESUB(L)+LYING_SNOW(L)/TIMESTEP TS1(L) = TS1(L) - SNOWMELT(L)*LF*TIMESTEP/HCAP_SOIL*DZ_SOIL else TS1(L) = TM endif endif !---------------------------------------------------------------------- ! Update the lying snow. !---------------------------------------------------------------------- LYING_SNOW(L) = LYING_SNOW(L) & + TIMESTEP*(SNOW(L)-ESUB(L)-SNOWMELT(L)) if (LYING_SNOW(L) .lt. 0.) then ! if negative snow, change extra sublimation to evaporation ESUB(L) = ESUB(L) + LYING_SNOW(L)/TIMESTEP E(L) = E(L) - LYING_SNOW(L)/TIMESTEP TS1(L) = TS1(L) + LF*LYING_SNOW(L)/(DZ_SOIL*HCAP_SOIL) LYING_SNOW(L) = 0. endif !---------------------------------------------------------------------- ! Update the soil moisture. !---------------------------------------------------------------------- RUNOFF(L) = KS*(M(L)/(1000.0*DZ_SOIL*THETA_SAT))**(2*B+3) M(L) = M(L) + TIMESTEP*(RAIN(L)+SNOWMELT(L)-E(L)-RUNOFF(L)) ! keep track of any negative moisture for conservation if (M(L) + MNEG(L) .lt. 0.) then MNEG(L) = MNEG(L) + M(L) M(L) = 0. else M(L) = MNEG(L) + M(L) MNEG(L) = 0. endif enddo # endif return end #endif