LCOV - code coverage report
Current view: top level - columnphysics - icepack_therm_0layer.F90 (source / functions) Hit Total Coverage
Test: 200704-015006:b1e41d9f12:3:base,travis,quick Lines: 49 85 57.65 %
Date: 2020-07-03 20:11:40 Functions: 1 1 100.00 %

          Line data    Source code
       1             : !=========================================================================
       2             : !
       3             : ! Update ice and snow internal temperatures
       4             : ! using zero-layer thermodynamics
       5             : !
       6             : ! authors: Alison McLaren, UK MetOffice
       7             : !          Elizabeth C. Hunke, LANL
       8             : !
       9             : ! 2012: Split from icepack_therm_vertical.F90
      10             : 
      11             :       module icepack_therm_0layer
      12             : 
      13             :       use icepack_kinds
      14             :       use icepack_parameters, only: c0, c1, p5, puny
      15             :       use icepack_parameters, only: kseaice, ksno
      16             :       use icepack_therm_bl99, only: surface_fluxes
      17             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      18             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      19             : 
      20             :       implicit none
      21             : 
      22             :       private
      23             :       public :: zerolayer_temperature
      24             : 
      25             : !=======================================================================
      26             : 
      27             :       contains
      28             : 
      29             : !=======================================================================
      30             : !
      31             : ! Compute new surface temperature using zero layer model of Semtner
      32             : ! (1976).
      33             : !
      34             : ! New temperatures are computed iteratively by solving a
      35             : ! surface flux balance equation (i.e. net surface flux from atmos
      36             : ! equals conductive flux from the top to the bottom surface).
      37             : !
      38             : ! author:  Alison McLaren, Met Office
      39             : !         (but largely taken from temperature_changes)
      40             : 
      41      303273 :       subroutine zerolayer_temperature(nilyr,    nslyr,    &
      42             :                                        rhoa,     flw,      &
      43             :                                        potT,     Qa,       &
      44             :                                        shcoef,   lhcoef,   &
      45             :                                        fswsfc,             &
      46             :                                        hilyr,    hslyr,    &
      47             :                                        Tsf,      Tbot,     &
      48             :                                        fsensn,   flatn,    &
      49             :                                        flwoutn,  fsurfn,   &
      50             :                                        fcondtopn,fcondbot  )
      51             : 
      52             :       integer (kind=int_kind), intent(in) :: &
      53             :          nilyr , & ! number of ice layers
      54             :          nslyr     ! number of snow layers
      55             : 
      56             :       real (kind=dbl_kind), intent(in) :: &
      57             :          rhoa        , & ! air density (kg/m^3)
      58             :          flw         , & ! incoming longwave radiation (W/m^2)
      59             :          potT        , & ! air potential temperature  (K)
      60             :          Qa          , & ! specific humidity (kg/kg)
      61             :          shcoef      , & ! transfer coefficient for sensible heat
      62             :          lhcoef      , & ! transfer coefficient for latent heat
      63             :          Tbot        , & ! ice bottom surface temperature (deg C)
      64             :          fswsfc          ! SW absorbed at ice/snow surface (W m-2)
      65             : 
      66             :       real (kind=dbl_kind), intent(in) :: &
      67             :          hilyr       , & ! ice layer thickness (m)
      68             :          hslyr           ! snow layer thickness (m)
      69             : 
      70             :       real (kind=dbl_kind), intent(inout):: &
      71             :          fsensn      , & ! surface downward sensible heat (W m-2)
      72             :          flatn       , & ! surface downward latent heat (W m-2)
      73             :          flwoutn     , & ! upward LW at surface (W m-2)
      74             :          fsurfn      , & ! net flux to top surface, excluding fcondtopn
      75             :          fcondtopn       ! downward cond flux at top surface (W m-2)
      76             : 
      77             :       real (kind=dbl_kind), intent(out):: &
      78             :          fcondbot        ! downward cond flux at bottom surface (W m-2)
      79             : 
      80             :       real (kind=dbl_kind), &
      81             :          intent(inout) :: &
      82             :          Tsf             ! ice/snow surface temperature, Tsfcn
      83             : 
      84             :       ! local variables
      85             : 
      86             :       logical (kind=log_kind), parameter :: &
      87             :          l_zerolayerchecks = .true.
      88             : 
      89             :       integer (kind=int_kind), parameter :: &
      90             :          nitermax = 50   ! max number of iterations in temperature solver
      91             : 
      92             :       real (kind=dbl_kind), parameter :: &
      93             :          Tsf_errmax = 5.e-4_dbl_kind ! max allowed error in Tsf
      94             :                                      ! recommend Tsf_errmax < 0.01 K
      95             : 
      96             :       integer (kind=int_kind) :: &
      97             :          niter           ! iteration counter in temperature solver
      98             : 
      99             :       real (kind=dbl_kind) :: &
     100      107286 :          Tsf_start   , & ! Tsf at start of iteration
     101      107286 :          dTsf        , & ! Tsf - Tsf_start
     102      107286 :          dfsurf_dT       ! derivative of fsurfn wrt Tsf
     103             : 
     104             :       real (kind=dbl_kind) :: &
     105      107286 :          dTsf_prev   , & ! dTsf from previous iteration
     106      107286 :          dfsens_dT   , & ! deriv of fsens wrt Tsf (W m-2 deg-1)
     107      107286 :          dflat_dT    , & ! deriv of flat wrt Tsf (W m-2 deg-1)
     108      107286 :          dflwout_dT      ! deriv of flwout wrt Tsf (W m-2 deg-1)
     109             : 
     110             :       real (kind=dbl_kind) :: &
     111      107286 :          kh          , & ! effective conductivity
     112      107286 :          diag        , & ! diagonal matrix elements
     113      107286 :          rhs             ! rhs of tri-diagonal matrix equation
     114             : 
     115             :       real (kind=dbl_kind) :: &
     116      107286 :          heff        , & ! effective ice thickness (m)
     117             :                          ! ( hice + hsno*kseaice/ksnow)
     118      107286 :          kratio      , & ! ratio of ice and snow conductivies
     119      107286 :          avg_Tsf         ! = 1. if Tsf averaged w/Tsf_start, else = 0.
     120             : 
     121             :       logical (kind=log_kind) :: &
     122             :          converged      ! = true when local solution has converged
     123             : 
     124             :       character(len=*),parameter :: subname='(zerolayer_temperature)'
     125             : 
     126             :       !-----------------------------------------------------------------
     127             :       ! Initialize
     128             :       !-----------------------------------------------------------------
     129             : 
     130      303273 :       fcondbot = c0
     131             : 
     132      303273 :       converged = .false.
     133             :       
     134      303273 :       dTsf_prev = c0
     135             : 
     136             :       !-----------------------------------------------------------------
     137             :       ! Solve for new temperatures.
     138             :       ! Iterate until temperatures converge with minimal temperature
     139             :       ! change.
     140             :       !-----------------------------------------------------------------
     141             : 
     142    15466923 :       do niter = 1, nitermax
     143             : 
     144    15466923 :          if (.not. converged) then
     145             :             
     146             :       !-----------------------------------------------------------------
     147             :       ! Update radiative and turbulent fluxes and their derivatives
     148             :       ! with respect to Tsf.
     149             :       !-----------------------------------------------------------------
     150             :             
     151             :             call surface_fluxes (Tsf,        fswsfc,            &
     152             :                                  rhoa,       flw,               &
     153             :                                  potT,       Qa,                &
     154             :                                  shcoef,     lhcoef,            &
     155             :                                  flwoutn,    fsensn,            &
     156             :                                  flatn,      fsurfn,            &
     157             :                                  dflwout_dT, dfsens_dT,         &
     158      417295 :                                  dflat_dT,   dfsurf_dT)
     159      417295 :             if (icepack_warnings_aborted(subname)) return
     160             : 
     161             :       !-----------------------------------------------------------------
     162             :       ! Compute effective ice thickness (includes snow) and thermal 
     163             :       ! conductivity 
     164             :       !-----------------------------------------------------------------
     165             : 
     166      417295 :             kratio = kseaice/ksno
     167             :    
     168      417295 :             heff = hilyr + kratio * hslyr
     169      417295 :             kh = kseaice / heff
     170             : 
     171             :       !-----------------------------------------------------------------
     172             :       ! Compute conductive flux at top surface, fcondtopn.
     173             :       ! If fsurfn < fcondtopn and Tsf = 0, then reset Tsf to slightly less
     174             :       !  than zero (but not less than -puny).
     175             :       !-----------------------------------------------------------------
     176             : 
     177      417295 :             fcondtopn = kh * (Tsf - Tbot)
     178             :             
     179      417295 :             if (fsurfn < fcondtopn) &
     180      378582 :                  Tsf = min (Tsf, -puny)
     181             :             
     182             :       !-----------------------------------------------------------------
     183             :       ! Save surface temperature at start of iteration
     184             :       !-----------------------------------------------------------------
     185             : 
     186      417295 :             Tsf_start = Tsf
     187             : 
     188             :       !-----------------------------------------------------------------
     189             :       ! Solve surface balance equation to obtain the new temperatures.
     190             :       !-----------------------------------------------------------------
     191             : 
     192      417295 :             diag = dfsurf_dT - kh
     193             :             rhs  = dfsurf_dT*Tsf - fsurfn   &
     194      417295 :                  - kh*Tbot
     195      417295 :             Tsf = rhs / diag
     196             : 
     197             :       !-----------------------------------------------------------------
     198             :       ! Determine whether the computation has converged to an acceptable
     199             :       ! solution.  Four conditions must be satisfied:
     200             :       !
     201             :       !    (1) Tsf <= 0 C.
     202             :       !    (2) Tsf is not oscillating; i.e., if both dTsf(niter) and
     203             :       !        dTsf(niter-1) have magnitudes greater than puny, then
     204             :       !        dTsf(niter)/dTsf(niter-1) cannot be a negative number
     205             :       !        with magnitude greater than 0.5.  
     206             :       !    (3) abs(dTsf) < Tsf_errmax
     207             :       !    (4) If Tsf = 0 C, then the downward turbulent/radiative 
     208             :       !        flux, fsurfn, must be greater than or equal to the downward
     209             :       !        conductive flux, fcondtopn.
     210             :       !-----------------------------------------------------------------
     211             : 
     212             :       !-----------------------------------------------------------------
     213             :       ! Initialize convergence flag (true until proven false), dTsf,
     214             :       !  and temperature-averaging coefficients.
     215             :       ! Average only if test 1 or 2 fails.
     216             :       ! Initialize energy.
     217             :       !-----------------------------------------------------------------
     218             : 
     219      417295 :             converged = .true.
     220      417295 :             dTsf = Tsf - Tsf_start
     221      417295 :             avg_Tsf = c0
     222             : 
     223             :       !-----------------------------------------------------------------
     224             :       ! Condition 1: check for Tsf > 0
     225             :       ! If Tsf > 0, set Tsf = 0 and leave converged=.true.
     226             :       !-----------------------------------------------------------------
     227             : 
     228      417295 :             if (Tsf > puny) then
     229           0 :                Tsf = c0
     230           0 :                dTsf = -Tsf_start
     231             : 
     232             :       !-----------------------------------------------------------------
     233             :       ! Condition 2: check for oscillating Tsf
     234             :       ! If oscillating, average all temps to increase rate of convergence.
     235             :       ! It is possible that this may never occur.
     236             :       !-----------------------------------------------------------------
     237             : 
     238             :             elseif (niter > 1 &                ! condition (2)
     239             :               .and. Tsf_start <= -puny &
     240             :               .and. abs(dTsf) > puny &
     241             :               .and. abs(dTsf_prev) > puny &
     242      417295 :               .and. -dTsf/(dTsf_prev+puny*puny) > p5) then
     243             : 
     244           0 :                avg_Tsf  = c1  ! average with starting temp  
     245           0 :                dTsf = p5 * dTsf
     246           0 :                converged = .false.
     247             :             endif
     248             : 
     249             :       !-----------------------------------------------------------------
     250             :       ! If condition 2 failed, average new surface temperature with
     251             :       !  starting value.
     252             :       !-----------------------------------------------------------------
     253             :             Tsf = Tsf &
     254      417295 :                 + avg_Tsf * p5 * (Tsf_start - Tsf)
     255             : 
     256             :       !-----------------------------------------------------------------
     257             :       ! Condition 3: check for large change in Tsf
     258             :       !-----------------------------------------------------------------
     259             : 
     260      417295 :             if (abs(dTsf) > Tsf_errmax) then
     261      114022 :                converged = .false.
     262             :             endif
     263             : 
     264             :       !-----------------------------------------------------------------
     265             :       ! Condition 4: check for fsurfn < fcondtopn with Tsf > 0
     266             :       !-----------------------------------------------------------------
     267             : 
     268      417295 :             fsurfn = fsurfn + dTsf*dfsurf_dT
     269      417295 :             fcondtopn = kh * (Tsf-Tbot)
     270             : 
     271      417295 :             if (Tsf > -puny .and. fsurfn < fcondtopn) then
     272           0 :                converged = .false.
     273             :             endif
     274             : 
     275      417295 :             fcondbot = fcondtopn
     276             : 
     277      417295 :             dTsf_prev = dTsf
     278             : 
     279             :          endif ! converged
     280             : 
     281             :       enddo                     ! temperature iteration niter
     282             : 
     283             :       !-----------------------------------------------------------------
     284             :       ! Check for convergence failures.
     285             :       !-----------------------------------------------------------------
     286      303273 :       if (.not.converged) then
     287           0 :          write(warnstr,*) subname, 'Thermo iteration does not converge,'
     288           0 :          call icepack_warnings_add(warnstr)
     289           0 :          write(warnstr,*) subname, 'Ice thickness:',  hilyr*nilyr
     290           0 :          call icepack_warnings_add(warnstr)
     291           0 :          write(warnstr,*) subname, 'Snow thickness:', hslyr*nslyr
     292           0 :          call icepack_warnings_add(warnstr)
     293           0 :          write(warnstr,*) subname, 'dTsf, Tsf_errmax:',dTsf_prev, &
     294           0 :                           Tsf_errmax
     295           0 :          call icepack_warnings_add(warnstr)
     296           0 :          write(warnstr,*) subname, 'Tsf:', Tsf
     297           0 :          call icepack_warnings_add(warnstr)
     298           0 :          write(warnstr,*) subname, 'fsurfn:', fsurfn
     299           0 :          call icepack_warnings_add(warnstr)
     300           0 :          write(warnstr,*) subname, 'fcondtopn, fcondbot', &
     301           0 :                           fcondtopn, fcondbot
     302           0 :          call icepack_warnings_add(warnstr)
     303           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     304           0 :          call icepack_warnings_add(subname//" zerolayer_temperature: Thermo iteration does not converge" ) 
     305           0 :          return
     306             :       endif
     307             : 
     308             :       !-----------------------------------------------------------------
     309             :       ! Check that if Tsfc < 0, then fcondtopn = fsurfn
     310             :       !-----------------------------------------------------------------
     311             : 
     312             :       if (l_zerolayerchecks) then
     313      303273 :          if (Tsf < c0 .and. & 
     314             :               abs(fcondtopn-fsurfn) > puny) then
     315             : 
     316           0 :             write(warnstr,*) subname, 'fcondtopn does not equal fsurfn,'
     317           0 :             call icepack_warnings_add(warnstr)
     318           0 :             write(warnstr,*) subname, 'Tsf=',Tsf
     319           0 :             call icepack_warnings_add(warnstr)
     320           0 :             write(warnstr,*) subname, 'fcondtopn=',fcondtopn
     321           0 :             call icepack_warnings_add(warnstr)
     322           0 :             write(warnstr,*) subname, 'fsurfn=',fsurfn
     323           0 :             call icepack_warnings_add(warnstr)
     324           0 :             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     325           0 :             call icepack_warnings_add(subname//" zerolayer_temperature: fcondtopn /= fsurfn" ) 
     326           0 :             return
     327             :          endif
     328             :       endif                     ! l_zerolayerchecks
     329             : 
     330             :       ! update fluxes that depend on Tsf
     331      303273 :       flwoutn = flwoutn + dTsf_prev * dflwout_dT
     332      303273 :       fsensn  = fsensn  + dTsf_prev * dfsens_dT
     333      303273 :       flatn   = flatn   + dTsf_prev * dflat_dT
     334             : 
     335             :       end subroutine zerolayer_temperature
     336             : 
     337             : !=======================================================================
     338             : 
     339             :       end module icepack_therm_0layer
     340             : 
     341             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd