fix more icepack debug errors

CICE-Consortium · Feb 6, 2018 · a971864 · a971864
1 parent 24755dd
commit a971864
Show file tree

Hide file tree

Showing 8 changed files with 47 additions and 34 deletions.
diff --git a/columnphysics/icepack_algae.F90 b/columnphysics/icepack_algae.F90
@@ -67,6 +67,9 @@ module icepack_algae
       private 
       public :: zbio, sklbio
 
+      real (kind=dbl_kind), parameter :: & 
+         exp_argmax = c10    ! maximum argument of exponential
+
 !=======================================================================
 
       contains
@@ -903,8 +906,8 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
          Sink_top,      & ! For cons: (+ or -) remaining bottom flux into ice(mmol/m^2/s)
          ocean_b,       & ! ocean_bio
          sum_react,     &
-         rtau_ret,      & ! retention frequency (s^-1)
-         rtau_rel     , & ! release frequency   (s^-1)
+         exp_ret,       & ! exp dt/retention frequency
+         exp_rel,       & ! exp dt/release frequency
          atm_add_cons , & ! zbgc_snow+zbgc_atm (mmol/m^3*m)
          dust_Fe      , & ! contribution of dust surface flux to dFe (umol/m*3*m)
          source       , & ! mmol/m^2 surface input from snow/atmosphere
@@ -933,7 +936,8 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
                           ! when hin > hbri:  just used in sw calculation
 
       real (kind=dbl_kind):: &
-         bio_tmp          ! temporary variable
+         bio_tmp, &       ! temporary variable
+         exp_min          ! temporary exp var
 
       real (kind=dbl_kind):: &
          Sat_conc   , & ! adsorbing saturation concentration  (mmols/m^3)
@@ -1058,21 +1062,27 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
       !-----------------------------------------------------------------
       !   time constants for mobile/stationary phase changes
       !-----------------------------------------------------------------
-
+
+         exp_rel(m) = c0
+         exp_ret(m) = c0
+         if (tau_ret(m) > c0) then
+            exp_min = min(dt/tau_ret(m),exp_argmax)
+            exp_ret(m) = exp(-exp_min)
+         endif
+         if (tau_rel(m) > c0) then
+            exp_min = min(dt/tau_rel(m),exp_argmax)
+            exp_rel(m) = exp(-exp_min)
+         endif
          if (m .ne. nlt_bgc_N(1)) then  
             if (hin_old  > hin) then  !melting
-               rtau_rel(m) = c1/tau_rel(m)
-               rtau_ret(m) = c0
+               exp_ret(m) = c1
             else                              !not melting
-               rtau_ret(m) = c1/tau_ret(m)
-               rtau_rel(m) = c0
+               exp_rel(m) = c1
             endif  
          elseif (tr_bgc_N .and. hin_old > hin + algal_vel*dt) then
-               rtau_rel(m) = c1/tau_rel(m)
-               rtau_ret(m) = c0
+               exp_ret(m) = c1
          elseif (tr_bgc_N) then
-               rtau_ret(m) = c1/tau_ret(m)
-               rtau_rel(m) = c0
+               exp_rel(m) = c1
          endif
 
          ocean_b(m) = ocean_bio(m)
@@ -1152,8 +1162,8 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
             do k = 1,nblyr+1
                initcons_mobile(k) = in_init_cons(k,mm)*trcrn(nt_zbgc_frac+mm-1)
                initcons_stationary(k) = mobile(mm)*(in_init_cons(k,mm)-initcons_mobile(k))
-               dmobile(k) = mobile(mm)*(initcons_mobile(k)*(exp(-dt*rtau_ret( mm))-c1) + &
-                                 initcons_stationary(k)*(c1-exp(-dt*rtau_rel(mm))))
+               dmobile(k) = mobile(mm)*(initcons_mobile(k)*(exp_ret(mm)-c1) + &
+                                    initcons_stationary(k)*(c1-exp_rel(mm)))
                initcons_mobile(k) = max(c0,initcons_mobile(k) + dmobile(k))
                initcons_stationary(k) = max(c0,initcons_stationary(k) - dmobile(k))
                if (initcons_stationary(k)/hbri_old > Sat_conc) then
@@ -1270,7 +1280,7 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
                 call icepack_warnings_add(warnstr)
                 write(warnstr,*) subname, 'in_init_cons(:,mm):',in_init_cons(:,mm)
                 call icepack_warnings_add(warnstr)
-                write(warnstr,*) subname, 'rtau_ret( mm),rtau_rel( mm)',rtau_ret( mm),rtau_rel( mm)
+                write(warnstr,*) subname, 'exp_ret( mm),exp_rel( mm)',exp_ret( mm),exp_rel( mm)
                 call icepack_warnings_add(warnstr)
                 write(warnstr,*) subname,'darcyV,dhtop,dhbot'
                 call icepack_warnings_add(warnstr)
@@ -1370,7 +1380,7 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
                 call icepack_warnings_add(warnstr)
                 write(warnstr,*) subname,  react(k,m),iphin_N(k),biomat_brine(k,m)
                 call icepack_warnings_add(warnstr)
-                write(warnstr,*) subname, 'rtau_ret( m),rtau_ret( m)',rtau_ret( m),rtau_ret( m)
+                write(warnstr,*) subname, 'exp_ret( m),exp_ret( m)',exp_ret( m),exp_ret( m)
                 call icepack_warnings_add(warnstr)
                 call icepack_warnings_add(subname//'negative bgc')
                 call icepack_warnings_setabort(.true.,__FILE__,__LINE__)

diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90
@@ -161,7 +161,7 @@ subroutine atmo_boundary_layer (sfctype,            &
          psimhu, & ! unstable part of psimh
          psixhu    ! unstable part of psimx
 
-      character(len=*),parameter :: subname='(atm_boundary_layer)'
+      character(len=*),parameter :: subname='(atmo_boundary_layer)'
 
       !------------------------------------------------------------
       ! Define functions
@@ -389,8 +389,7 @@ subroutine atmo_boundary_const (sfctype,  calc_strair, &
                                       Tsf,      potT,     &
                                       Qa,                 &
                                       delt,     delq,     &
-                                      lhcoef,   shcoef,   &
-                                      Cdn_atm)  
+                                      lhcoef,   shcoef    )
 
       character (len=3), intent(in) :: &
          sfctype      ! ice or ocean
@@ -407,9 +406,6 @@ subroutine atmo_boundary_const (sfctype,  calc_strair, &
          wind     , & ! wind speed (m/s)
          rhoa         ! air density (kg/m^3)
 
-      real (kind=dbl_kind), intent(in) :: &
-         Cdn_atm      ! neutral drag coefficient
-
       real (kind=dbl_kind), intent(inout):: &
          strx     , & ! x surface stress (N)
          stry         ! y surface stress (N)
@@ -429,7 +425,7 @@ subroutine atmo_boundary_const (sfctype,  calc_strair, &
          tau, &  ! stress at zlvl
          Lheat   ! Lvap or Lsub, depending on surface type
 
-      character(len=*),parameter :: subname='(atm_boundary_const)'
+      character(len=*),parameter :: subname='(atmo_boundary_const)'
 
       !------------------------------------------------------------
       ! Initialize
@@ -901,6 +897,8 @@ subroutine icepack_atm_boundary(sfctype,                    &
          worku = uvel
       endif
 
+      Cdn_atm_ratio_n = c1
+
                if (trim(atmbndy) == 'constant') then
                   call atmo_boundary_const (sfctype,  calc_strair, &
                                             uatm,     vatm,     &
@@ -909,8 +907,7 @@ subroutine icepack_atm_boundary(sfctype,                    &
                                             Tsf,      potT,     &
                                             Qa,                 &
                                             delt,     delq,     &
-                                            lhcoef,   shcoef,   &
-                                            Cdn_atm)
+                                            lhcoef,   shcoef    )
                   if (icepack_warnings_aborted(subname)) return
                else ! default
                   call atmo_boundary_layer (sfctype,                 &

diff --git a/columnphysics/icepack_brine.F90 b/columnphysics/icepack_brine.F90
@@ -8,7 +8,7 @@
       module icepack_brine
 
       use icepack_kinds
-      use icepack_parameters, only: p01, p001, p5, c0, c1, c2, c1p5, puny
+      use icepack_parameters, only: p01, p001, p5, c0, c1, c2, c1p5, c10, puny
       use icepack_parameters, only: gravit, rhoi, rhow, rhos, depressT
       use icepack_parameters, only: salt_loss, min_salin, rhosi
       use icepack_parameters, only: dts_b, l_sk
@@ -45,6 +45,9 @@ module icepack_brine
          b1      = 1000.0_dbl_kind, & ! (kg/m^3)  
          b2      = 0.8_dbl_kind       ! (kg/m^3/ppt)
 
+      real (kind=dbl_kind), parameter :: & 
+         exp_argmax = c10    ! maximum argument of exponential
+
 !=======================================================================
 
       contains
@@ -511,6 +514,7 @@ subroutine update_hbrine (meltb,      meltt,       &
          darcy_coeff, & ! magnitude of the Darcy velocity/hbrocn (1/s)
          hbrocn_new , & ! hbrocn after flushing
          dhflood    , & ! surface flooding by ocean
+         exp_min    , & ! temporary exp value
          dhrunoff       ! direct runoff to ocean
 
       real (kind=dbl_kind), parameter :: &
@@ -540,7 +544,8 @@ subroutine update_hbrine (meltb,      meltt,       &
                bphi_min   = bphin  
                dhrunoff  = -dhS_top*aice0
                hbrocn    = max(c0,hbrocn - dhrunoff)
-               hbrocn_new = hbrocn*exp(-darcy_coeff/bphi_min*dt)
+               exp_min = min(darcy_coeff/bphi_min*dt,exp_argmax)
+               hbrocn_new = hbrocn*exp(-exp_min)
                hbr = max(hbrmin, h_ocn + hbrocn_new)
                hbrocn_new = hbr-h_ocn
                darcy_V = -SIGN((hbrocn-hbrocn_new)/dt*bphi_min, hbrocn)

diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
@@ -678,7 +678,7 @@ subroutine init_vertical_profile(nilyr,    nslyr,    &
          zSin            ! internal ice layer salinities
 
       real (kind=dbl_kind), dimension (:), &
-         intent(out) :: &
+         intent(inout) :: &
          zqsn        , & ! snow enthalpy
          zTsn            ! snow temperature
 

diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90
@@ -135,7 +135,7 @@ subroutine init_restart
       use icedrv_init_column, only: init_hbrine, init_bgc
       use icedrv_restart, only: restartfile
       use icedrv_restart_shared, only: restart
-      use icedrv_restart_column, only: read_restart_bgc
+      use icedrv_restart_bgc, only: read_restart_bgc
       use icedrv_state ! almost everything
 
       integer(kind=int_kind) :: &

diff --git a/configuration/driver/icedrv_RunMod.F90 b/configuration/driver/icedrv_RunMod.F90
@@ -94,7 +94,7 @@ subroutine ice_step
       use icedrv_flux, only: scale_factor, init_history_therm, init_history_bgc, &
           daidtt, daidtd, dvidtt, dvidtd, dagedtt, dagedtd, init_history_dyn
       use icedrv_restart, only: dumpfile, final_restart
-      use icedrv_restart_column, only: write_restart_bgc
+      use icedrv_restart_bgc, only: write_restart_bgc
       use icedrv_state, only: trcrn
       use icedrv_step, only: prep_radiation, step_therm1, step_therm2, &
           update_state, step_dyn_ridge, step_radiation, &

diff --git a/...guration/driver/icedrv_restart_column.F90 → configuration/driver/icedrv_restart_bgc.F90 b/...guration/driver/icedrv_restart_column.F90 → configuration/driver/icedrv_restart_bgc.F90
@@ -4,7 +4,7 @@
 !
 ! author: Elizabeth C. Hunke, LANL
 !
-      module icedrv_restart_column
+      module icedrv_restart_bgc
 
       use icedrv_kinds
       use icedrv_constants
@@ -20,7 +20,8 @@ module icedrv_restart_column
       implicit none
 
       private
-      public ::  write_restart_bgc,       read_restart_bgc
+      public ::  write_restart_bgc, &
+                 read_restart_bgc
 
 !=======================================================================
 
@@ -667,6 +668,6 @@ end subroutine read_restart_bgc
 
 !=======================================================================
 
-      end module icedrv_restart_column
+      end module icedrv_restart_bgc
 
 !=======================================================================
diff --git a/doc/source/developer_guide/dg_driver.rst b/doc/source/developer_guide/dg_driver.rst
@@ -29,7 +29,7 @@ The icepack driver consists of the following files
 |        **icedrv_init.F90**        general initialization routines
 |        **icedrv_init_column.F90** initialization routines specific to the column physics
 |        **icedrv_restart.F90**     driver for reading/writing restart files
-|        **icedrv_restart_column.F90**  restart routines specific to the column physics
+|        **icedrv_restart_bgc.F90**  restart routines specific to the column physics
 |        **icedrv_restart_shared.F90**  code shared by all restart options
 |        **icedrv_state.F90**       essential arrays to describe the state of the ice
 |        **icedrv_step.F90**        routines for time stepping the major code components