Merge pull request #78 from gustavo-marques/update_melt_potential

Improve melt potential calculation

config_src/coupled_driver/ocean_model_MOM.F90

@@ -240,6 +240,12 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn)
   ! Local variables
   real :: Rho0        ! The Boussinesq ocean density, in kg m-3.
   real :: G_Earth     ! The gravitational acceleration in m s-2.
+  real :: HFrz        !< If HFrz > 0 (m), melt potential will be computed.
+                      !! The actual depth over which melt potential is computed will
+                      !! min(HFrz, OBLD), where OBLD is the boundary layer depth.
+                      !! If HFrz <= 0 (default), melt potential will not be computed.
+  logical :: use_melt_pot!< If true, allocate melt_potential array
+
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mdl = "ocean_model_init"  ! This module's name.
@@ -342,8 +348,20 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn)
 
   !   Consider using a run-time flag to determine whether to do the diagnostic
   ! vertical integrals, since the related 3-d sums are not negligible in cost.
-  call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, &
-                              do_integrals=.true., gas_fields_ocn=gas_fields_ocn)
+  call get_param(param_file, mdl, "HFREEZE", HFrz, &
+                 "If HFREEZE > 0, melt potential will be computed. The actual depth \n"//&
+                 "over which melt potential is computed will be min(HFREEZE, OBLD), \n"//&
+                 "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), \n"//&
+                 "melt potential will not be computed.", units="m", default=-1.0, do_not_log=.true.)
+
+  if (HFrz .gt. 0.0) then
+    use_melt_pot=.true.
+  else
+    use_melt_pot=.false.
+  endif
+
+  call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, do_integrals=.true., &
+                              gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot)
 
   call surface_forcing_init(Time_in, OS%grid, param_file, OS%diag, &
                             OS%forcing_CSp, OS%restore_salinity, OS%restore_temp)

config_src/mct_driver/MOM_ocean_model.F90

@@ -243,9 +243,15 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
 ! Because of the way that indicies and domains are handled, Ocean_sfc must have
 ! been used in a previous call to initialize_ocean_type.
 
-  real :: Rho0        !< The Boussinesq ocean density, in kg m-3.
-  real :: G_Earth     !< The gravitational acceleration in m s-2.
-                      !! This include declares and sets the variable "version".
+  real :: Rho0           !< The Boussinesq ocean density, in kg m-3.
+  real :: G_Earth        !< The gravitational acceleration in m s-2.
+                         !! This include declares and sets the variable "version".
+  real :: HFrz           !< If HFrz > 0 (m), melt potential will be computed.
+                         !! The actual depth over which melt potential is computed will
+                         !! min(HFrz, OBLD), where OBLD is the boundary layer depth.
+                         !! If HFrz <= 0 (default), melt potential will not be computed.
+  logical :: use_melt_pot!< If true, allocate melt_potential array
+
 #include "version_variable.h"
   character(len=40)  :: mdl = "ocean_model_init"  !< This module's name.
   character(len=48)  :: stagger
@@ -338,8 +344,21 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
 
   !   Consider using a run-time flag to determine whether to do the diagnostic
   ! vertical integrals, since the related 3-d sums are not negligible in cost.
+
+  call get_param(param_file, mdl, "HFREEZE", HFrz, &
+                 "If HFREEZE > 0, melt potential will be computed. The actual depth \n"//&
+                 "over which melt potential is computed will be min(HFREEZE, OBLD), \n"//&
+                 "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), \n"//&
+                 "melt potential will not be computed.", units="m", default=-1.0, do_not_log=.true.)
+
+  if (HFrz .gt. 0.0) then
+    use_melt_pot=.true.
+  else
+    use_melt_pot=.false.
+  endif
+
   call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, do_integrals=.true., &
-                              gas_fields_ocn=gas_fields_ocn, use_meltpot=.true.)
+                              gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot)
 
   call surface_forcing_init(Time_in, OS%grid, param_file, OS%diag, &
                             OS%forcing_CSp, OS%restore_salinity, OS%restore_temp)

config_src/mct_driver/ocn_comp_mct.F90

@@ -310,7 +310,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename )
   glb%grid => glb%ocn_state%grid
 
   ! Allocate IOB data type (needs to be called after glb%grid is set)
-  write(6,*)'DEBUG: isc,iec,jsc,jec= ',glb%grid%isc, glb%grid%iec, glb%grid%jsc, glb%grid%jec
+  !write(6,*)'DEBUG: isc,iec,jsc,jec= ',glb%grid%isc, glb%grid%iec, glb%grid%jsc, glb%grid%jec
   call IOB_allocate(ice_ocean_boundary, glb%grid%isc, glb%grid%iec, glb%grid%jsc, glb%grid%jec)
 
   call t_stopf('MOM_init')

src/core/MOM.F90

@@ -275,6 +275,10 @@ module MOM
                                 !! average surface tracer properties (in meter) when
                                 !! bulk mixed layer is not used, or a negative value
                                 !! if a bulk mixed layer is being used.
+  real :: HFrz                  !< If HFrz > 0, melt potential will be computed.
+                                !! The actual depth over which melt potential is computed will
+                                !! min(HFrz, OBLD), where OBLD is the boundary layer depth.
+                                !! If HFrz <= 0 (default), melt potential will not be computed.
   real :: Hmix_UV               !< Depth scale over which to average surface flow to
                                 !! feedback to the coupler/driver (m) when
                                 !! bulk mixed layer is not used, or a negative value
@@ -1757,6 +1761,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
                  "SSU, SSV. A non-positive value indicates no averaging.", &
                  units="m", default=0.)
   endif
+  call get_param(param_file, "MOM", "HFREEZE", CS%HFrz, &
+                 "If HFREEZE > 0, melt potential will be computed. The actual depth \n"//&
+                 "over which melt potential is computed will be min(HFREEZE, OBLD), \n"//&
+                 "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), \n"//&
+                 "melt potential will not be computed.", units="m", default=-1.0)
   call get_param(param_file, "MOM", "MIN_Z_DIAG_INTERVAL", Z_diag_int, &
                  "The minimum amount of time in seconds between \n"//&
                  "calculations of depth-space diagnostics. Making this \n"//&
@@ -2673,6 +2682,7 @@ subroutine extract_surface_state(CS, sfc_state)
   real :: dh                          !< thickness of a layer within mixed layer (meter)
   real :: mass                        !< mass per unit area of a layer (kg/m2)
   real :: T_freeze                    !< freezing temperature (oC)
+  real :: delT(SZI_(CS%G))            !< T-T_freeze (oC)
   logical :: use_temperature   !< If true, temp and saln used as state variables.
   integer :: i, j, k, is, ie, js, je, nz, numberOfErrors
   integer :: isd, ied, jsd, jed
@@ -2831,21 +2841,41 @@ subroutine extract_surface_state(CS, sfc_state)
     endif
   endif  ! (CS%Hmix >= 0.0)
 
+
   if (allocated(sfc_state%melt_potential)) then
-    !$OMP parallel do default(shared)
-    do j=js,je ; do i=is,ie
-      ! set melt_potential to zero to avoid passing values set previously
-      if (G%mask2dT(i,j)>0.) then
-        ! calculate freezing pot. temp. @ surface
-        call calculate_TFreeze(sfc_state%SSS(i,j), 0.0, T_freeze, CS%tv%eqn_of_state)
-        ! time accumulated melt_potential, in J/m^2
-        sfc_state%melt_potential(i,j) = sfc_state%melt_potential(i,j) +  (CS%tv%C_p * CS%GV%Rho0 * &
-                                        (sfc_state%SST(i,j) - T_freeze) * CS%Hmix)
-      else
-          sfc_state%melt_potential(i,j) = 0.0
-      endif! G%mask2dT
-    enddo ; enddo
-  endif
+  !$OMP parallel do default(shared)
+    do j=js,je
+      do i=is,ie
+        depth(i) = 0.0
+        delT(i) = 0.0
+      enddo
+
+      do k=1,nz ; do i=is,ie
+        depth_ml = min(CS%HFrz,CS%visc%MLD(i,j))
+        if (depth(i) + h(i,j,k)*GV%H_to_m < depth_ml) then
+          dh = h(i,j,k)*GV%H_to_m
+        elseif (depth(i) < depth_ml) then
+          dh = depth_ml - depth(i)
+        else
+          dh = 0.0
+        endif
+
+        ! p=0 OK, HFrz ~ 10 to 20m
+        call calculate_TFreeze(CS%tv%S(i,j,k), 0.0, T_freeze, CS%tv%eqn_of_state)
+        depth(i) = depth(i) + dh
+        delT(i) =  delT(i) + dh * (CS%tv%T(i,j,k) - T_freeze)
+      enddo ; enddo
+
+      do i=is,ie
+       if (G%mask2dT(i,j)>0.) then
+         ! time accumulated melt_potential, in J/m^2
+         sfc_state%melt_potential(i,j) = sfc_state%melt_potential(i,j) +  (CS%tv%C_p * CS%GV%Rho0 * delT(i))
+       else
+         sfc_state%melt_potential(i,j) = 0.0
+       endif! G%mask2dT
+      enddo
+    enddo ! end of j loop
+  endif   ! melt_potential
 
   if (allocated(sfc_state%salt_deficit) .and. associated(CS%tv%salt_deficit)) then
     !$OMP parallel do default(shared)

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -615,6 +615,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
     if (associated(Hml)) then
       call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G)
       call pass_var(Hml, G%domain, halo=1)
+      ! If visc%MLD exists, copy KPP's BLD into it
+      if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:)
     endif
 
     call cpu_clock_end(id_clock_kpp)
@@ -1532,6 +1534,8 @@ subroutine legacy_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_en
     if (associated(Hml)) then
       call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G)
       call pass_var(Hml, G%domain, halo=1)
+      ! If visc%MLD exists, copy KPP's BLD into it
+      if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:)
     endif
 
     if (.not. CS%KPPisPassive) then

src/parameterizations/vertical/MOM_set_viscosity.F90

@@ -1723,6 +1723,7 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS)
   logical :: adiabatic, useKPP, useEPBL
   logical :: use_CVMix_shear, MLE_use_PBL_MLD, use_CVMix_conv
   integer :: isd, ied, jsd, jed, nz
+  real :: hfreeze !< If hfreeze > 0 (m), melt potential will be computed.
   character(len=40)  :: mdl = "MOM_set_visc"  ! This module's name.
   isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke
 
@@ -1779,15 +1780,24 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS)
                 "Vertical turbulent viscosity at interfaces due to slow processes", &
                 "m2 s-1", z_grid='i')
 
-  ! visc%MLD is used to communicate the state of the (e)PBL to the rest of the model
+  ! visc%MLD is used to communicate the state of the (e)PBL or KPP to the rest of the model
   call get_param(param_file, mdl, "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, &
                  default=.false., do_not_log=.true.)
+  ! visc%MLD needs to be allocated when melt potential is computed (HFREEZE>0)
+  call get_param(param_file, mdl, "HFREEZE", hfreeze, &
+                 default=-1.0, do_not_log=.true.)
+
   if (MLE_use_PBL_MLD) then
     call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
     call register_restart_field(visc%MLD, "MLD", .false., restart_CS, &
                   "Instantaneous active mixing layer depth", "m")
   endif
 
+  if (hfreeze >= 0.0 .and. .not.MLE_use_PBL_MLD) then
+    call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
+  endif
+
+
 end subroutine set_visc_register_restarts
 
 !> Initializes the MOM_set_visc control structure