From 37833eb7ed044d021f8529888f02c9bc2ca386bb Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 24 May 2016 16:39:52 -0400 Subject: [PATCH 01/46] Changes in the sponge layer Sponge layer is now initialized in the same way as the ocean interior. Most sponge layer parameters are now specified in the MOM_input file. --- src/user/ISOMIP_initialization.F90 | 167 +++++++++++++++++++++-------- 1 file changed, 124 insertions(+), 43 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index ba3ec38a81..d022ab64cc 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -181,6 +181,8 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) e0(k) = -G%max_depth * ( 0.5 * ( GV%Rlay(k-1) + GV%Rlay(k) ) - rho_sur ) / rho_range e0(k) = min( 0., e0(k) ) ! Bound by surface e0(k) = max( -G%max_depth, e0(k) ) ! Bound by possible deepest point in model + write(*,*)'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k) + enddo e0(nz+1) = -G%max_depth @@ -244,7 +246,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, ! real :: T_ref, T_Light, T_Dense, S_ref, S_Light, S_Dense, a1, frac_dense, k_frac, res_rat, k_light real :: z ! vertical position in z space character(len=40) :: verticalCoordinate, density_profile -real :: rho_light, delta_rho, z_tmp, rho_tmp + real :: rho_tmp is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -260,16 +262,13 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, do_not_log=.true.) call calculate_density(t_sur,s_sur,0.0,rho_sur,eqn_of_state) - write (*,*)'Density in the surface layer:', rho_sur + !write (*,*)'Density in the surface layer:', rho_sur call calculate_density(t_bot,s_bot,0.0,rho_bot,eqn_of_state) - write (*,*)'Density in the bottom layer::', rho_bot - - call get_param(param_file, mod, "LIGHTEST_DENSITY", rho_light) - call get_param(param_file, mod, "DENSITY_RANGE", delta_rho) + !write (*,*)'Density in the bottom layer::', rho_bot S_range = s_sur - s_bot T_range = t_sur - t_bot - write(*,*)'s_bot, s_sur, S_range, t_bot, t_sur, T_range', s_bot, s_sur, S_range, t_bot, t_sur, T_range + !write(*,*)'s_bot, s_sur, S_range, t_bot, t_sur, T_range', s_bot, s_sur, S_range, t_bot, t_sur, T_range select case ( coordinateMode(verticalCoordinate) ) @@ -288,7 +287,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, i=G%iec; j=G%jec do k = 1,nz call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,eqn_of_state) - write(0,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) + write(*,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) enddo case default @@ -324,50 +323,117 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) real :: RHO(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for RHO real :: h(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for thickness real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate, in s-1. - real :: S_ref, T_ref; ! Reference salinity and temerature within - ! surface layer - real :: S_range, T_range; ! Range of salinities and temperatures over the - ! vertical - + real :: TNUDG ! Nudging time scale, days + real :: S_sur, T_sur; ! Surface salinity and temerature in sponge + real :: S_bot, T_bot; ! Bottom salinity and temerature in sponge + real :: t_ref, s_ref ! reference T and S + real :: rho_sur, rho_bot, rho_range, t_range, s_range real :: e0(SZK_(G)) ! The resting interface heights, in m, usually ! ! negative because it is positive upward. ! real :: eta1D(SZK_(G)+1) ! Interface height relative to the sea surface ! ! positive upward, in m. - real :: min_depth, dummy1, z, Bmax - real :: damp, rho_dummy + real :: min_depth, dummy1, z, delta_h + real :: damp, rho_dummy, min_thickness, rho_tmp, xi0 + character(len=40) :: verticalCoordinate + character(len=40) :: mod = "ISOMIP_initialize_sponges" ! This subroutine's name. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + call get_param(PF,mod,"MIN_THICKNESS",min_thickness,'Minimum layer thickness',units='m',default=1.e-3) + + call get_param(PF,mod,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & + default=DEFAULT_COORDINATE_MODE) + + call get_param(PF, mod, "TNUDG", TNUDG, 'Nudging time scale for sponge layers (days)', default=0.0) + + call get_param(PF, mod, "T_REF", t_ref, 'Reference temperature', default=10.0,& + do_not_log=.true.) + + call get_param(PF, mod, "S_REF", s_ref, 'Reference salinity', default=35.0,& + do_not_log=.true.) + + call get_param(PF, mod, "S_SUR_SPONGE", s_sur, 'Surface salinity in sponge layer.', default=s_ref) + + call get_param(PF, mod, "S_BOT_SPONGE", s_bot, 'Bottom salinity in sponge layer.', default=s_ref) + + call get_param(PF, mod, "T_SUR_SPONGE", t_sur, 'Surface temperature in sponge layer.', default=t_ref) + + call get_param(PF, mod, "T_BOT_SPONGE", t_bot, 'Bottom temperature in sponge layer.', default=t_ref) + T(:,:,:) = 0.0 ; S(:,:,:) = 0.0 ; Idamp(:,:) = 0.0; RHO(:,:,:) = 0.0 - S_ref = 33.8; S_range = 0.90 - T_ref = -1.9; T_range = 2.9 - Bmax = 720.0 + S_range = s_sur - s_bot + T_range = t_sur - t_bot ! Set up sponges for ISOMIP configuration call get_param(PF, mod, "MINIMUM_DEPTH", min_depth, & "The minimum depth of the ocean.", units="m", default=0.0) -! GMM: set thickness, I am using same config. as the initial thickness for now - do k=1,nz - e0(k) = -G%max_depth * real(k-1) / real(nz) - enddo - - do j=js,je ; do i=is,ie - eta1D(nz+1) = -1.0*G%bathyT(i,j) - do k=nz,1,-1 - eta1D(k) = e0(k) - if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_z)) then - eta1D(k) = eta1D(k+1) + GV%Angstrom_z - h(i,j,k) = GV%Angstrom_z - else - h(i,j,k) = eta1D(k) - eta1D(k+1) - endif - enddo - enddo; enddo + + select case ( coordinateMode(verticalCoordinate) ) + + case ( REGRIDDING_LAYER, REGRIDDING_RHO ) + ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT + call calculate_density(t_sur,s_sur,0.0,rho_sur,tv%eqn_of_state) + write (*,*)'Surface density in sponge:', rho_sur + call calculate_density(t_bot,s_bot,0.0,rho_bot,tv%eqn_of_state) + write (*,*)'Bottom density in sponge:', rho_bot + rho_range = rho_bot - rho_sur + write (*,*)'Density range in sponge:', rho_range + + ! Construct notional interface positions + e0(1) = 0. + do K=2,nz + e0(k) = -G%max_depth * ( 0.5 * ( GV%Rlay(k-1) + GV%Rlay(k) ) - rho_sur ) / rho_range + e0(k) = min( 0., e0(k) ) ! Bound by surface + e0(k) = max( -G%max_depth, e0(k) ) ! Bound by possible deepest point in model + !write(*,*)'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k) + + enddo + e0(nz+1) = -G%max_depth + + ! Calculate thicknesses + do j=js,je ; do i=is,ie + eta1D(nz+1) = -1.0*G%bathyT(i,j) + do k=nz,1,-1 + eta1D(k) = e0(k) + if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_z)) then + eta1D(k) = eta1D(k+1) + GV%Angstrom_z + h(i,j,k) = GV%Angstrom_z + else + h(i,j,k) = eta1D(k) - eta1D(k+1) + endif + enddo + enddo ; enddo + + case ( REGRIDDING_ZSTAR ) ! Initial thicknesses for z coordinates + do j=js,je ; do i=is,ie + eta1D(nz+1) = -1.0*G%bathyT(i,j) + do k=nz,1,-1 + eta1D(k) = -G%max_depth * real(k-1) / real(nz) + if (eta1D(k) < (eta1D(k+1) + min_thickness)) then + eta1D(k) = eta1D(k+1) + min_thickness + h(i,j,k) = min_thickness + else + h(i,j,k) = eta1D(k) - eta1D(k+1) + endif + enddo + enddo ; enddo + + case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates + do j=js,je ; do i=is,ie + delta_h = G%bathyT(i,j) / dfloat(nz) + h(i,j,:) = delta_h + end do ; end do + + case default + call MOM_error(FATAL,"ISOMIP_initialize_sponges: "// & + "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE") + + end select ! Here the inverse damping time, in s-1, is set. Set Idamp to 0 ! ! wherever there is no sponge, and the subroutines that are called ! @@ -379,8 +445,7 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) ! 1 / day dummy1=(G%geoLonT(i,j)-790.0)/(800.0-790.0) -! damp = 1.0/10.0 * max(0.0,dummy1) - damp = 1.0/0.04166 * max(0.0,dummy1) ! one hour + damp = 1.0/TNUDG * max(0.0,dummy1) else ; damp=0.0 endif @@ -405,16 +470,32 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) call initialize_ALE_sponge(Idamp,h, nz, G, PF, CSp) ! setup temp and salt at the sponge zone + select case ( coordinateMode(verticalCoordinate) ) + + case ( REGRIDDING_SIGMA, REGRIDDING_ZSTAR, REGRIDDING_RHO ) + S_range = S_range / G%max_depth ! Convert S_range into dS/dz + T_range = T_range / G%max_depth ! Convert T_range into dT/dz do j=js,je ; do i=is,ie + xi0 = -G%bathyT(i,j); do k = nz,1,-1 -! z = (G%bathyT(i,j)/(nz-1))* (k -1) - z = (G%max_depth/(nz-1))* (k -1) - S(i,j,k) = S_REF + (S_RANGE*z/Bmax) - T(i,j,k) = T_REF + (T_RANGE*z/Bmax) -! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_dummy,tv%eqn_of_state) -! RHO(i,j,k) = rho_dummy + xi0 = xi0 + 0.5 * h(i,j,k) ! Depth in middle of layer + S(i,j,k) = S_sur + S_range * xi0 + T(i,j,k) = T_sur + T_range * xi0 + xi0 = xi0 + 0.5 * h(i,j,k) ! Depth at top of layer enddo enddo ; enddo + i=G%iec; j=G%jec + do k = 1,nz + call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state) + write(*,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) + enddo + + case default + call MOM_error(FATAL,"isomip_initialize_sponge: "// & + "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE") + + end select + ! Now register all of the fields which are damped in the sponge. ! ! By default, momentum is advected vertically within the sponge, but ! From ce0b2cf6685020f310feebd9eab62eef60d81dab Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 24 May 2016 20:11:09 -0400 Subject: [PATCH 02/46] Added a passive tracer to the ISOMIP test case Created a module that sets up and use (one for now) a passive tracers that is injected in the sponge layer. --- src/tracer/ISOMIP_tracer.F90 | 412 +++++++++++++++++++++++++ src/tracer/MOM_tracer_flow_control.F90 | 20 ++ 2 files changed, 432 insertions(+) create mode 100644 src/tracer/ISOMIP_tracer.F90 diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 new file mode 100644 index 0000000000..48d25d3bec --- /dev/null +++ b/src/tracer/ISOMIP_tracer.F90 @@ -0,0 +1,412 @@ +!> This module contains the routines used to set up and use a set of (one for now) +!! dynamically passive tracers. For now, just one passive tracer is injected in +!! the sponge layer. +!! Set up and use passive tracers requires the following: +!! (1) register_ISOMIP_tracer +!! (2) + +!********+*********+*********+*********+*********+*********+*********+** +!* * +!* By Robert Hallberg, 2002 * +!* Adapted to the ISOMIP test case by Gustavo Marques, May 2016 * !* * +!********+*********+*********+*********+*********+*********+*********+** + +module ISOMIP_tracer + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr +use MOM_diag_mediator, only : diag_ctrl +use MOM_diag_to_Z, only : register_Z_tracer, diag_to_Z_CS +use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_forcing_type, only : forcing +use MOM_grid, only : ocean_grid_type +use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_restart, only : register_restart_field, MOM_restart_CS +use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS +use MOM_time_manager, only : time_type, get_time +use MOM_tracer_registry, only : register_tracer, tracer_registry_type +use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values +use MOM_tracer_registry, only : tracer_vertdiff +use MOM_variables, only : surface, ocean_OBC_type +use MOM_verticalGrid, only : verticalGrid_type + +use coupler_util, only : set_coupler_values, ind_csurf +use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux + +implicit none ; private + +#include + +!< Publicly available functions +public register_ISOMIP_tracer, initialize_ISOMIP_tracer +public ISOMIP_tracer_column_physics, ISOMIP_tracer_surface_state, ISOMIP_tracer_end + +!< ntr is the number of tracers in this module. +integer, parameter :: ntr = 1 + +type p3d + real, dimension(:,:,:), pointer :: p => NULL() +end type p3d + +!> tracer control structure +type, public :: ISOMIP_tracer_CS ; private + logical :: coupled_tracers = .false. !< These tracers are not offered to the + !< coupler. + character(len = 200) :: tracer_IC_file !< The full path to the IC file, or " " + !< to initialize internally. + type(time_type), pointer :: Time !< A pointer to the ocean model's clock. + type(tracer_registry_type), pointer :: tr_Reg => NULL() + real, pointer :: tr(:,:,:,:) => NULL() !< The array of tracers used in this + !< subroutine, in g m-3? + real, pointer :: tr_aux(:,:,:,:) => NULL() !< The masked tracer concentration + !< for output, in g m-3. + type(p3d), dimension(NTR) :: & + tr_adx, &!< Tracer zonal advective fluxes in g m-3 m3 s-1. + tr_ady, &!< Tracer meridional advective fluxes in g m-3 m3 s-1. + tr_dfx, &!< Tracer zonal diffusive fluxes in g m-3 m3 s-1. + tr_dfy !< Tracer meridional diffusive fluxes in g m-3 m3 s-1. + real :: land_val(NTR) = -1.0 !< The value of tr used where land is masked out. + logical :: mask_tracers !< If true, tracers are masked out in massless layers. + logical :: use_sponge + + integer, dimension(NTR) :: ind_tr !< Indices returned by aof_set_coupler_flux + !< if it is used and the surface tracer concentrations are to be + !< provided to the coupler. + + type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the + !< timing of diagnostic output. + integer, dimension(NTR) :: id_tracer = -1, id_tr_adx = -1, id_tr_ady = -1 + integer, dimension(NTR) :: id_tr_dfx = -1, id_tr_dfy = -1 + + type(vardesc) :: tr_desc(NTR) +end type ISOMIP_tracer_CS + +contains + + +!> This subroutine is used to register tracer fields +function register_ISOMIP_tracer(G, param_file, CS, diag, tr_Reg, & + restart_CS) + type(ocean_grid_type), intent(in) :: G ! NULL() + logical :: register_ISOMIP_tracer + integer :: isd, ied, jsd, jed, nz, m + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = G%ke + + if (associated(CS)) then + call MOM_error(WARNING, "ISOMIP_register_tracer called with an "// & + "associated control structure.") + return + endif + allocate(CS) + + CS%diag => diag + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mod, version, "") + call get_param(param_file, mod, "ISOMIP_TRACER_IC_FILE", CS%tracer_IC_file, & + "The name of a file from which to read the initial \n"//& + "conditions for the ISOMIP tracers, or blank to initialize \n"//& + "them internally.", default=" ") + if (len_trim(CS%tracer_IC_file) >= 1) then + call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + CS%tracer_IC_file = trim(inputdir)//trim(CS%tracer_IC_file) + call log_param(param_file, mod, "INPUTDIR/ISOMIP_TRACER_IC_FILE", & + CS%tracer_IC_file) + endif + call get_param(param_file, mod, "SPONGE", CS%use_sponge, & + "If true, sponges may be applied anywhere in the domain. \n"//& + "The exact location and properties of those sponges are \n"//& + "specified from MOM_initialization.F90.", default=.false.) + + allocate(CS%tr(isd:ied,jsd:jed,nz,NTR)) ; CS%tr(:,:,:,:) = 0.0 + if (CS%mask_tracers) then + allocate(CS%tr_aux(isd:ied,jsd:jed,nz,NTR)) ; CS%tr_aux(:,:,:,:) = 0.0 + endif + + do m=1,NTR + if (m < 10) then ; write(name,'("tr_D",I1.1)') m + else ; write(name,'("tr_D",I2.2)') m ; endif + write(longname,'("Concentration of ISOMIP Tracer ",I2.2)') m + CS%tr_desc(m) = var_desc(name, units="kg kg-1", longname=longname, caller=mod) + + ! This is needed to force the compiler not to do a copy in the registration + ! calls. Curses on the designers and implementers of Fortran90. + tr_ptr => CS%tr(:,:,:,m) + ! Register the tracer for the restart file. + call register_restart_field(tr_ptr, CS%tr_desc(m), .true., restart_CS) + ! Register the tracer for horizontal advection & diffusion. + call register_tracer(tr_ptr, CS%tr_desc(m), param_file, G, tr_Reg, & + tr_desc_ptr=CS%tr_desc(m)) + + ! Set coupled_tracers to be true (hard-coded above) to provide the surface + ! values to the coupler (if any). This is meta-code and its arguments will + ! currently (deliberately) give fatal errors if it is used. + if (CS%coupled_tracers) & + CS%ind_tr(m) = aof_set_coupler_flux(trim(name)//'_flux', & + flux_type=' ', implementation=' ', caller="register_ISOMIP_tracer") + enddo + + CS%tr_Reg => tr_Reg + register_ISOMIP_tracer = .true. +end function register_ISOMIP_tracer + +!> Initializes the NTR tracer fields in tr(:,:,:,:) +! and it sets up the tracer output. +subroutine initialize_ISOMIP_tracer(restart, day, G, GV, h, OBC, CS, ALE_sponge_CSp, & + diag_to_Z_CSp) + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + logical, intent(in) :: restart !< .true. if the fields have already been read from a restart file. + type(time_type), target, intent(in) :: day !< Time of the start of the run. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, where, and what open boundary conditions are used. This is not being used for now. + type(ISOMIP_tracer_CS), pointer :: CS !< The control structure returned by a previous call to ISOMIP_register_tracer. + type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< A pointer to the control structure for the sponges, if they are in use. Otherwise this may be unassociated. + type(diag_to_Z_CS), pointer :: diag_to_Z_CSp !< A pointer to the control structure for diagnostics in depth space. + + real, allocatable :: temp(:,:,:) + real, pointer, dimension(:,:,:) :: & + OBC_tr1_u => NULL(), & ! These arrays should be allocated and set to + OBC_tr1_v => NULL() ! specify the values of tracer 1 that should come + ! in through u- and v- points through the open + ! boundary conditions, in the same units as tr. + character(len=16) :: name ! A variable's name in a NetCDF file. + character(len=72) :: longname ! The long name of that variable. + character(len=48) :: units ! The dimensions of the variable. + character(len=48) :: flux_units ! The units for tracer fluxes, usually + ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1. + real, pointer :: tr_ptr(:,:,:) => NULL() + real :: PI ! 3.1415926... calculated as 4*atan(1) + real :: tr_y ! Initial zonally uniform tracer concentrations. + real :: dist2 ! The distance squared from a line, in m2. + real :: h_neglect ! A thickness that is so small it is usually lost + ! in roundoff and can be neglected, in m. + real :: e(SZK_(G)+1), e_top, e_bot, d_tr + integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m + integer :: IsdB, IedB, JsdB, JedB + + if (.not.associated(CS)) return + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + h_neglect = GV%H_subroundoff + + CS%Time => day + + if (.not.restart) then + if (len_trim(CS%tracer_IC_file) >= 1) then + ! Read the tracer concentrations from a netcdf file. + if (.not.file_exists(CS%tracer_IC_file, G%Domain)) & + call MOM_error(FATAL, "ISOMIP_initialize_tracer: Unable to open "// & + CS%tracer_IC_file) + do m=1,NTR + call query_vardesc(CS%tr_desc(m), name, caller="initialize_ISOMIP_tracer") + call read_data(CS%tracer_IC_file, trim(name), & + CS%tr(:,:,:,m), domain=G%Domain%mpp_domain) + enddo + else + do m=1,NTR + do k=1,nz ; do j=js,je ; do i=is,ie + CS%tr(i,j,k,m) = 1.0e-20 ! This could just as well be 0. + enddo ; enddo ; enddo + enddo + endif + endif ! restart + + if ( CS%use_sponge ) then +! If sponges are used, this example damps tracers in sponges in the +! northern half of the domain to 1 and tracers in the southern half +! to 0. For any tracers that are not damped in the sponge, the call +! to set_up_sponge_field can simply be omitted. + if (.not.associated(ALE_sponge_CSp)) & + call MOM_error(FATAL, "ISOMIP_initialize_tracer: "// & + "The pointer to ALEsponge_CSp must be associated if SPONGE is defined.") + + allocate(temp(G%isd:G%ied,G%jsd:G%jed,nz)) + + do j=js,je ; do i=is,ie + if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then + temp(i,j,:) = 1.0 + else + temp(i,j,:) = 0.0 + endif + enddo ; enddo + +! do m=1,NTR + do m=1,1 + ! This is needed to force the compiler not to do a copy in the sponge + ! calls. Curses on the designers and implementers of Fortran90. + tr_ptr => CS%tr(:,:,:,m) + call set_up_ALE_sponge_field(temp, G, tr_ptr, ALE_sponge_CSp) + enddo + deallocate(temp) + endif + + ! This needs to be changed if the units of tracer are changed above. + if (GV%Boussinesq) then ; flux_units = "kg kg-1 m3 s-1" + else ; flux_units = "kg s-1" ; endif + + do m=1,NTR + ! Register the tracer for the restart file. + call query_vardesc(CS%tr_desc(m), name, units=units, longname=longname, & + caller="initialize_ISOMIP_tracer") + CS%id_tracer(m) = register_diag_field("ocean_model", trim(name), CS%diag%axesTL, & + day, trim(longname) , trim(units)) + CS%id_tr_adx(m) = register_diag_field("ocean_model", trim(name)//"_adx", & + CS%diag%axesCuL, day, trim(longname)//" advective zonal flux" , & + trim(flux_units)) + CS%id_tr_ady(m) = register_diag_field("ocean_model", trim(name)//"_ady", & + CS%diag%axesCvL, day, trim(longname)//" advective meridional flux" , & + trim(flux_units)) + CS%id_tr_dfx(m) = register_diag_field("ocean_model", trim(name)//"_dfx", & + CS%diag%axesCuL, day, trim(longname)//" diffusive zonal flux" , & + trim(flux_units)) + CS%id_tr_dfy(m) = register_diag_field("ocean_model", trim(name)//"_dfy", & + CS%diag%axesCvL, day, trim(longname)//" diffusive zonal flux" , & + trim(flux_units)) + if (CS%id_tr_adx(m) > 0) call safe_alloc_ptr(CS%tr_adx(m)%p,IsdB,IedB,jsd,jed,nz) + if (CS%id_tr_ady(m) > 0) call safe_alloc_ptr(CS%tr_ady(m)%p,isd,ied,JsdB,JedB,nz) + if (CS%id_tr_dfx(m) > 0) call safe_alloc_ptr(CS%tr_dfx(m)%p,IsdB,IedB,jsd,jed,nz) + if (CS%id_tr_dfy(m) > 0) call safe_alloc_ptr(CS%tr_dfy(m)%p,isd,ied,JsdB,JedB,nz) + +! Register the tracer for horizontal advection & diffusion. + if ((CS%id_tr_adx(m) > 0) .or. (CS%id_tr_ady(m) > 0) .or. & + (CS%id_tr_dfx(m) > 0) .or. (CS%id_tr_dfy(m) > 0)) & + call add_tracer_diagnostics(name, CS%tr_Reg, CS%tr_adx(m)%p, & + CS%tr_ady(m)%p, CS%tr_dfx(m)%p, CS%tr_dfy(m)%p) + + call register_Z_tracer(CS%tr(:,:,:,m), trim(name), longname, units, & + day, G, diag_to_Z_CSp) + enddo + +end subroutine initialize_ISOMIP_tracer + +!> This subroutine applies diapycnal diffusion and any other column +! tracer physics or chemistry to the tracers from this file. +! This is a simple example of a set of advected passive tracers. +subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS) + type(ocean_grid_type), intent(in) :: G + type(verticalGrid_type), intent(in) :: GV + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old, h_new, ea, eb + type(forcing), intent(in) :: fluxes + real, intent(in) :: dt + type(ISOMIP_tracer_CS), pointer :: CS + +! Arguments: h_old - Layer thickness before entrainment, in m or kg m-2. +! (in) h_new - Layer thickness after entrainment, in m or kg m-2. +! (in) ea - an array to which the amount of fluid entrained +! from the layer above during this call will be +! added, in m or kg m-2. +! (in) eb - an array to which the amount of fluid entrained +! from the layer below during this call will be +! added, in m or kg m-2. +! (in) fluxes - A structure containing pointers to any possible +! forcing fields. Unused fields have NULL ptrs. +! (in) dt - The amount of time covered by this call, in s. +! (in) G - The ocean's grid structure. +! (in) GV - The ocean's vertical grid structure. +! (in) CS - The control structure returned by a previous call to +! ISOMIP_register_tracer. +! +! The arguments to this subroutine are redundant in that +! h_new[k] = h_old[k] + ea[k] - eb[k-1] + eb[k] - ea[k+1] + + real :: b1(SZI_(G)) ! b1 and c1 are variables used by the + real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver. + integer :: i, j, k, is, ie, js, je, nz, m + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + + if (.not.associated(CS)) return + + do m=1,NTR + call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV) + enddo + + if (CS%mask_tracers) then + do m = 1,NTR ; if (CS%id_tracer(m) > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + if (h_new(i,j,k) < 1.1*GV%Angstrom) then + CS%tr_aux(i,j,k,m) = CS%land_val(m) + else + CS%tr_aux(i,j,k,m) = CS%tr(i,j,k,m) + endif + enddo ; enddo ; enddo + endif ; enddo + endif + + do m=1,NTR + if (CS%mask_tracers) then + if (CS%id_tracer(m)>0) & + call post_data(CS%id_tracer(m),CS%tr_aux(:,:,:,m),CS%diag) + else + if (CS%id_tracer(m)>0) & + call post_data(CS%id_tracer(m),CS%tr(:,:,:,m),CS%diag) + endif + if (CS%id_tr_adx(m)>0) & + call post_data(CS%id_tr_adx(m),CS%tr_adx(m)%p(:,:,:),CS%diag) + if (CS%id_tr_ady(m)>0) & + call post_data(CS%id_tr_ady(m),CS%tr_ady(m)%p(:,:,:),CS%diag) + if (CS%id_tr_dfx(m)>0) & + call post_data(CS%id_tr_dfx(m),CS%tr_dfx(m)%p(:,:,:),CS%diag) + if (CS%id_tr_dfy(m)>0) & + call post_data(CS%id_tr_dfy(m),CS%tr_dfy(m)%p(:,:,:),CS%diag) + enddo + +end subroutine ISOMIP_tracer_column_physics + +!> This particular tracer package does not report anything back to the coupler. +! The code that is here is just a rough guide for packages that would. +subroutine ISOMIP_tracer_surface_state(state, h, G, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(surface), intent(inout) :: state !< A structure containing fields that describe the surface state of the ocean. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(ISOMIP_tracer_CS), pointer :: CS !< The control structure returned by a previous call to ISOMIP_register_tracer. + integer :: i, j, m, is, ie, js, je, nz + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + + if (.not.associated(CS)) return + + if (CS%coupled_tracers) then + do m=1,ntr + ! This call loads the surface vlues into the appropriate array in the + ! coupler-type structure. + call set_coupler_values(CS%tr(:,:,1,1), state%tr_fields, CS%ind_tr(m), & + ind_csurf, is, ie, js, je) + enddo + endif + +end subroutine ISOMIP_tracer_surface_state + +subroutine ISOMIP_tracer_end(CS) + type(ISOMIP_tracer_CS), pointer :: CS + integer :: m + + if (associated(CS)) then + if (associated(CS%tr)) deallocate(CS%tr) + if (associated(CS%tr_aux)) deallocate(CS%tr_aux) + do m=1,NTR + if (associated(CS%tr_adx(m)%p)) deallocate(CS%tr_adx(m)%p) + if (associated(CS%tr_ady(m)%p)) deallocate(CS%tr_ady(m)%p) + if (associated(CS%tr_dfx(m)%p)) deallocate(CS%tr_dfx(m)%p) + if (associated(CS%tr_dfy(m)%p)) deallocate(CS%tr_dfy(m)%p) + enddo + + deallocate(CS) + endif +end subroutine ISOMIP_tracer_end + +end module ISOMIP_tracer diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index b5df8c7151..b46e5aa083 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -51,6 +51,9 @@ module MOM_tracer_flow_control use DOME_tracer, only : register_DOME_tracer, initialize_DOME_tracer use DOME_tracer, only : DOME_tracer_column_physics, DOME_tracer_surface_state use DOME_tracer, only : DOME_tracer_end, DOME_tracer_CS +use ISOMIP_tracer, only : register_ISOMIP_tracer, initialize_ISOMIP_tracer +use ISOMIP_tracer, only : ISOMIP_tracer_column_physics, ISOMIP_tracer_surface_state +use ISOMIP_tracer, only : ISOMIP_tracer_end, ISOMIP_tracer_CS use ideal_age_example, only : register_ideal_age_tracer, initialize_ideal_age_tracer use ideal_age_example, only : ideal_age_tracer_column_physics, ideal_age_tracer_surface_state use ideal_age_example, only : ideal_age_stock, ideal_age_example_end, ideal_age_tracer_CS @@ -79,6 +82,7 @@ module MOM_tracer_flow_control type, public :: tracer_flow_control_CS ; private logical :: use_USER_tracer_example = .false. logical :: use_DOME_tracer = .false. + logical :: use_ISOMIP_tracer = .false. logical :: use_ideal_age = .false. logical :: use_oil = .false. logical :: use_advection_test_tracer = .false. @@ -86,6 +90,7 @@ module MOM_tracer_flow_control logical :: use_MOM_generic_tracer = .false. type(USER_tracer_example_CS), pointer :: USER_tracer_example_CSp => NULL() type(DOME_tracer_CS), pointer :: DOME_tracer_CSp => NULL() + type(ISOMIP_tracer_CS), pointer :: ISOMIP_tracer_CSp => NULL() type(ideal_age_tracer_CS), pointer :: ideal_age_tracer_CSp => NULL() type(oil_tracer_CS), pointer :: oil_tracer_CSp => NULL() type(advection_test_tracer_CS), pointer :: advection_test_tracer_CSp => NULL() @@ -137,6 +142,9 @@ subroutine call_tracer_register(G, param_file, CS, diag, tr_Reg, restart_CS) call get_param(param_file, mod, "USE_DOME_TRACER", CS%use_DOME_tracer, & "If true, use the DOME_tracer tracer package.", & default=.false.) + call get_param(param_file, mod, "USE_ISOMIP_TRACER", CS%use_ISOMIP_tracer, & + "If true, use the ISOMIP_tracer tracer package.", & + default=.false.) call get_param(param_file, mod, "USE_IDEAL_AGE_TRACER", CS%use_ideal_age, & "If true, use the ideal_age_example tracer package.", & default=.false.) @@ -169,6 +177,9 @@ subroutine call_tracer_register(G, param_file, CS, diag, tr_Reg, restart_CS) if (CS%use_DOME_tracer) CS%use_DOME_tracer = & register_DOME_tracer(G, param_file, CS%DOME_tracer_CSp, & diag, tr_Reg, restart_CS) + if (CS%use_ISOMIP_tracer) CS%use_ISOMIP_tracer = & + register_ISOMIP_tracer(G, param_file, CS%ISOMIP_tracer_CSp, & + diag, tr_Reg, restart_CS) if (CS%use_ideal_age) CS%use_ideal_age = & register_ideal_age_tracer(G, param_file, CS%ideal_age_tracer_CSp, & diag, tr_Reg, restart_CS) @@ -230,6 +241,9 @@ subroutine tracer_flow_control_init(restart, day, G, GV, h, param_file, OBC, & if (CS%use_DOME_tracer) & call initialize_DOME_tracer(restart, day, G, GV, h, OBC, CS%DOME_tracer_CSp, & sponge_CSp, diag_to_Z_CSp) + if (CS%use_ISOMIP_tracer) & + call initialize_ISOMIP_tracer(restart, day, G, GV, h, OBC, CS%ISOMIP_tracer_CSp, & + ALE_sponge_CSp, diag_to_Z_CSp) if (CS%use_ideal_age) & call initialize_ideal_age_tracer(restart, day, G, GV, h, OBC, CS%ideal_age_tracer_CSp, & sponge_CSp, diag_to_Z_CSp) @@ -346,6 +360,9 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, dt, G, GV, tv, o if (CS%use_DOME_tracer) & call DOME_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & G, GV, CS%DOME_tracer_CSp) + if (CS%use_ISOMIP_tracer) & + call ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, CS%ISOMIP_tracer_CSp) if (CS%use_ideal_age) & call ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & G, GV, CS%ideal_age_tracer_CSp) @@ -532,6 +549,8 @@ subroutine call_tracer_surface_state(state, h, G, CS) call USER_tracer_surface_state(state, h, G, CS%USER_tracer_example_CSp) if (CS%use_DOME_tracer) & call DOME_tracer_surface_state(state, h, G, CS%DOME_tracer_CSp) + if (CS%use_ISOMIP_tracer) & + call ISOMIP_tracer_surface_state(state, h, G, CS%ISOMIP_tracer_CSp) if (CS%use_ideal_age) & call ideal_age_tracer_surface_state(state, h, G, CS%ideal_age_tracer_CSp) if (CS%use_oil) & @@ -553,6 +572,7 @@ subroutine tracer_flow_control_end(CS) if (CS%use_USER_tracer_example) & call USER_tracer_example_end(CS%USER_tracer_example_CSp) if (CS%use_DOME_tracer) call DOME_tracer_end(CS%DOME_tracer_CSp) + if (CS%use_ISOMIP_tracer) call ISOMIP_tracer_end(CS%ISOMIP_tracer_CSp) if (CS%use_ideal_age) call ideal_age_example_end(CS%ideal_age_tracer_CSp) if (CS%use_oil) call oil_tracer_end(CS%oil_tracer_CSp) if (CS%use_advection_test_tracer) call advection_test_tracer_end(CS%advection_test_tracer_CSp) From af08bc0e0dd73403595e63c3b94e27eec96064f6 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 25 May 2016 17:28:18 -0400 Subject: [PATCH 03/46] Added a passive tracer to the ice shelf melt water I am now adding a second passive tracer (m=2) based on the ice shelf melt water field. To do so I first included the melt field in the structure that contains pointers to the boundary forcing (defined in MOM_forcing_type.F90) so this field can be accessed from ISOMIP_tracer.F90. This field is passed to forcing structure in MOM_ice_shelf.F90. Finally, the tracer is set in ISOMIP_tracer_column_physics (ISOMIP_tracer.F90). PS: this needs to be improved. It is still not working the way it should. --- src/core/MOM_forcing_type.F90 | 26 ++++++++++++++++++-------- src/ice_shelf/MOM_ice_shelf.F90 | 6 ++++++ src/tracer/ISOMIP_tracer.F90 | 13 ++++++++++++- 3 files changed, 36 insertions(+), 9 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index e96888250a..e2f9998b6a 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -121,14 +121,15 @@ module MOM_forcing_type ! land ice-shelf related inputs real, pointer, dimension(:,:) :: & - ustar_shelf => NULL(), & !< friction velocity under ice-shelves (m/s) - !! as computed by the ocean at the previous time step. - frac_shelf_h => NULL(), & !< Fractional ice shelf coverage of h-, u-, and v- - frac_shelf_u => NULL(), & !< cells, nondimensional from 0 to 1. These are only - frac_shelf_v => NULL(), & !< associated if ice shelves are enabled, and are - !! exactly 0 away from shelves or on land. - rigidity_ice_u => NULL(),& !< Depth-integrated lateral viscosity of ice - rigidity_ice_v => NULL() !< shelves or sea ice at u- or v-points (m3/s) + ustar_shelf => NULL(), & !< friction velocity under ice-shelves (m/s) + !! as computed by the ocean at the previous time step. + frac_shelf_h => NULL(), & !< Fractional ice shelf coverage of h-, u-, and v- + frac_shelf_u => NULL(), & !< cells, nondimensional from 0 to 1. These are only + frac_shelf_v => NULL(), & !< associated if ice shelves are enabled, and are + !! exactly 0 away from shelves or on land. + iceshelf_melt => NULL(), & !< ice shelf melt rate (positive) or freezing (negative) ( m/year ) + rigidity_ice_u => NULL(),& !< Depth-integrated lateral viscosity of ice + rigidity_ice_v => NULL() !< shelves or sea ice at u- or v-points (m3/s) ! Scalars set by surface forcing modules real :: vPrecGlobalAdj !< adjustment to restoring vprec to zero out global net ( kg/(m^2 s) ) @@ -1615,6 +1616,13 @@ subroutine forcing_accumulate(flux_tmp, fluxes, dt, G, wt2) fluxes%ustar_shelf(i,j) = flux_tmp%ustar_shelf(i,j) enddo ; enddo endif + + if (associated(fluxes%iceshelf_melt) .and. associated(flux_tmp%iceshelf_melt)) then + do i=isd,ied ; do j=jsd,jed + fluxes%iceshelf_melt(i,j) = flux_tmp%iceshelf_melt(i,j) + enddo ; enddo + endif + if (associated(fluxes%frac_shelf_h) .and. associated(flux_tmp%frac_shelf_h)) then do i=isd,ied ; do j=jsd,jed fluxes%frac_shelf_h(i,j) = flux_tmp%frac_shelf_h(i,j) @@ -2198,6 +2206,7 @@ subroutine allocate_forcing_type(G, fluxes, stress, ustar, water, heat, shelf, p call myAlloc(fluxes%frac_shelf_u,IsdB,IedB,jsd,jed, shelf) call myAlloc(fluxes%frac_shelf_v,isd,ied,JsdB,JedB, shelf) call myAlloc(fluxes%ustar_shelf,isd,ied,jsd,jed, shelf) + call myAlloc(fluxes%iceshelf_melt,isd,ied,jsd,jed, shelf) call myAlloc(fluxes%rigidity_ice_u,IsdB,IedB,jsd,jed, shelf) call myAlloc(fluxes%rigidity_ice_v,isd,ied,JsdB,JedB, shelf) @@ -2263,6 +2272,7 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%TKE_tidal)) deallocate(fluxes%TKE_tidal) if (associated(fluxes%ustar_tidal)) deallocate(fluxes%ustar_tidal) if (associated(fluxes%ustar_shelf)) deallocate(fluxes%ustar_shelf) + if (associated(fluxes%iceshelf_melt)) deallocate(fluxes%iceshelf_melt) if (associated(fluxes%frac_shelf_h)) deallocate(fluxes%frac_shelf_h) if (associated(fluxes%frac_shelf_u)) deallocate(fluxes%frac_shelf_u) if (associated(fluxes%frac_shelf_v)) deallocate(fluxes%frac_shelf_v) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 79dfb09a40..44109551f9 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -682,6 +682,9 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) enddo ! i-loop enddo ! j-loop + ! melt in m/year + fluxes%iceshelf_melt = CS%lprec * (86400.0*365.0/CS%density_ice) + if (CS%DEBUG) then call hchksum (CS%h_shelf, "melt rate", G, haloshift=0) endif @@ -1382,6 +1385,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti allocate( fluxes%frac_shelf_u(Isdq:Iedq,jsd:jed) ) ; fluxes%frac_shelf_u(:,:) = 0.0 allocate( fluxes%frac_shelf_v(isd:ied,Jsdq:Jedq) ) ; fluxes%frac_shelf_v(:,:) = 0.0 allocate( fluxes%ustar_shelf(isd:ied,jsd:jed) ) ; fluxes%ustar_shelf(:,:) = 0.0 + allocate( fluxes%iceshelf_melt(isd:ied,jsd:jed) ) ; fluxes%iceshelf_melt(:,:) = 0.0 if (.not.associated(fluxes%p_surf)) then allocate( fluxes%p_surf(isd:ied,jsd:jed) ) ; fluxes%p_surf(:,:) = 0.0 endif @@ -1445,6 +1449,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti if (.not. solo_ice_sheet) then vd = var_desc("ustar_shelf","m s-1","Friction velocity under ice shelves",z_grid='1') call register_restart_field(fluxes%ustar_shelf, vd, .true., CS%restart_CSp) + vd = var_desc("iceshelf_melt","m year-1","Ice Shelf Melt Rate",z_grid='1') + call register_restart_field(fluxes%iceshelf_melt, vd, .true., CS%restart_CSp) endif CS%restart_output_dir = dirs%restart_output_dir diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index 48d25d3bec..3fdf4374ef 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -44,7 +44,7 @@ module ISOMIP_tracer public ISOMIP_tracer_column_physics, ISOMIP_tracer_surface_state, ISOMIP_tracer_end !< ntr is the number of tracers in this module. -integer, parameter :: ntr = 1 +integer, parameter :: ntr = 2 type p3d real, dimension(:,:,:), pointer :: p => NULL() @@ -336,6 +336,17 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV) enddo + ! dye melt water (m=2) + ! converting melt (m/year) to (m/s) + do m=2,NTR + do j=js,je ; do i=is,ie + if (fluxes%iceshelf_melt(i,j) > 0.0) then + CS%tr(i,j,1,m) = (dt*fluxes%iceshelf_melt(i,j)/(365.0 * 86400.0) & + + h_old(i,j,1)*CS%tr(i,j,1,m))/h_new(i,j,1) + endif + enddo ; enddo + enddo + if (CS%mask_tracers) then do m = 1,NTR ; if (CS%id_tracer(m) > 0) then do k=1,nz ; do j=js,je ; do i=is,ie From 33e088fce06e5dc22accae08108e9f35faa9030d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 27 May 2016 12:35:18 -0400 Subject: [PATCH 04/46] Changed names of parameters Changed the name of parameters that are exclusively to the ISOMIP test case. --- src/ice_shelf/MOM_ice_shelf.F90 | 1 - src/tracer/ISOMIP_tracer.F90 | 12 ++-- src/user/ISOMIP_initialization.F90 | 93 ++++++++++++------------------ 3 files changed, 43 insertions(+), 63 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 44109551f9..e375497ae3 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -678,7 +678,6 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) else !not shelf CS%t_flux(i,j) = 0.0 endif - enddo ! i-loop enddo ! j-loop diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index 3fdf4374ef..152980b04b 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -332,21 +332,23 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G if (.not.associated(CS)) return - do m=1,NTR - call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV) - enddo ! dye melt water (m=2) ! converting melt (m/year) to (m/s) do m=2,NTR do j=js,je ; do i=is,ie if (fluxes%iceshelf_melt(i,j) > 0.0) then - CS%tr(i,j,1,m) = (dt*fluxes%iceshelf_melt(i,j)/(365.0 * 86400.0) & - + h_old(i,j,1)*CS%tr(i,j,1,m))/h_new(i,j,1) +! CS%tr(i,j,1,m) = (dt*fluxes%iceshelf_melt(i,j)/(365.0 * 86400.0) & +! + h_old(i,j,1)*CS%tr(i,j,1,m))/h_new(i,j,1) + CS%tr(i,j,1,m) = 1.0 endif enddo ; enddo enddo + do m=1,NTR + call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV) + enddo + if (CS%mask_tracers) then do m = 1,NTR ; if (CS%id_tracer(m) > 0) then do k=1,nz ; do j=js,je ; do i=is,ie diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index d022ab64cc..7c8264f74b 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -129,18 +129,12 @@ end subroutine ISOMIP_initialize_topography !> Initialization of thicknesses subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) - type(ocean_grid_type), intent(in) :: G - type(verticalGrid_type), intent(in) :: GV - real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h - type(param_file_type), intent(in) :: param_file - type(thermo_var_ptrs), intent(in) :: tv - -! Arguments: h - The thickness that is being initialized. -! (in) G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) tv - A structure containing pointers to any available -! thermodynamic fields, including eq. of state + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The container for vertical grid data + real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is being initialized. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to parse for model parameter values. + type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available thermodynamic fields, including eq. of state. + real :: e0(SZK_(G)+1) ! The resting interface heights, in m, usually ! ! negative because it is positive upward. ! real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface ! @@ -162,18 +156,18 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) select case ( coordinateMode(verticalCoordinate) ) case ( REGRIDDING_LAYER, REGRIDDING_RHO ) ! Initial thicknesses for isopycnal coordinates - call get_param(param_file, mod, "T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9) - call get_param(param_file, mod, "S_SUF", s_sur, 'Salinity at the surface (interface)', default=33.8) - call get_param(param_file, mod, "T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9) - call get_param(param_file, mod, "S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55) + call get_param(param_file, mod, "ISOMIP_T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9) + call get_param(param_file, mod, "ISOMIP_S_SUR", s_sur, 'Salinity at the surface (interface)', default=33.8) + call get_param(param_file, mod, "ISOMIP_T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9) + call get_param(param_file, mod, "ISOMIP_S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55) ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT call calculate_density(t_sur,s_sur,0.0,rho_sur,tv%eqn_of_state) - write (*,*)'Surface density is:', rho_sur + !write (*,*)'Surface density is:', rho_sur call calculate_density(t_bot,s_bot,0.0,rho_bot,tv%eqn_of_state) - write (*,*)'Bottom density is:', rho_bot + !write (*,*)'Bottom density is:', rho_bot rho_range = rho_bot - rho_sur - write (*,*)'Density range is:', rho_range + !write (*,*)'Density range is:', rho_range ! Construct notional interface positions e0(1) = 0. @@ -181,7 +175,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) e0(k) = -G%max_depth * ( 0.5 * ( GV%Rlay(k-1) + GV%Rlay(k) ) - rho_sur ) / rho_range e0(k) = min( 0., e0(k) ) ! Bound by surface e0(k) = max( -G%max_depth, e0(k) ) ! Bound by possible deepest point in model - write(*,*)'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k) + !write(*,*)'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k) enddo e0(nz+1) = -G%max_depth @@ -243,7 +237,6 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, integer :: i, j, k, is, ie, js, je, nz real :: x, ds, dt, rho_sur, rho_bot real :: xi0, xi1, dxi, r, S_sur, T_sur, S_bot, T_bot, S_range, T_range -! real :: T_ref, T_Light, T_Dense, S_ref, S_Light, S_Dense, a1, frac_dense, k_frac, res_rat, k_light real :: z ! vertical position in z space character(len=40) :: verticalCoordinate, density_profile real :: rho_tmp @@ -252,13 +245,13 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, call get_param(param_file,mod,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & default=DEFAULT_COORDINATE_MODE) - call get_param(param_file, mod, "T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9,& + call get_param(param_file, mod, "ISOMIP_T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9,& do_not_log=.true.) - call get_param(param_file, mod, "S_SUF", s_sur, 'Salinity at the surface (interface)', default=33.8,& + call get_param(param_file, mod, "ISOMIP_S_SUR", s_sur, 'Salinity at the surface (interface)', default=33.8,& do_not_log=.true.) - call get_param(param_file, mod, "T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9,& + call get_param(param_file, mod, "ISOMIP_T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9,& do_not_log=.true.) - call get_param(param_file, mod, "S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55,& + call get_param(param_file, mod, "ISOMIP_S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55,& do_not_log=.true.) call calculate_density(t_sur,s_sur,0.0,rho_sur,eqn_of_state) @@ -268,7 +261,6 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, S_range = s_sur - s_bot T_range = t_sur - t_bot - !write(*,*)'s_bot, s_sur, S_range, t_bot, t_sur, T_range', s_bot, s_sur, S_range, t_bot, t_sur, T_range select case ( coordinateMode(verticalCoordinate) ) @@ -287,7 +279,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, i=G%iec; j=G%jec do k = 1,nz call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,eqn_of_state) - write(*,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) + !write(*,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) enddo case default @@ -302,22 +294,12 @@ end subroutine ISOMIP_initialize_temperature_salinity ! the values towards which the interface heights and an arbitrary ! ! number of tracers should be restored within each sponge. subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) - type(ocean_grid_type), intent(in) :: G - type(verticalGrid_type), intent(in) :: GV - type(thermo_var_ptrs), intent(in) :: tv - type(param_file_type), intent(in) :: PF - type(ALE_sponge_CS), pointer :: CSp - -! Arguments: G - The ocean's grid structure. -! (in) tv - A structure containing pointers to any available -! thermodynamic fields, including potential temperature and -! salinity or mixed layer density. Absent fields have NULL ptrs. -! (in) PF - A structure indicating the open file to parse for -! model parameter values. -! (in/out) CSp - A pointer that is set to point to the control structure -! for this module - -! real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1) ! A temporary array for eta. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The container for vertical grid data + type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available thermodynamic fields, including potential temperature and salinity or mixed layer density. Absent fields have NULL ptrs. + type(param_file_type), intent(in) :: PF !< A structure indicating the open file to parse for model parameter values. + type(ALE_sponge_CS), pointer :: CSp !< A pointer that is set to point to the control structure for this module. + real :: T(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for temp real :: S(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for salt real :: RHO(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for RHO @@ -348,21 +330,21 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) call get_param(PF,mod,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & default=DEFAULT_COORDINATE_MODE) - call get_param(PF, mod, "TNUDG", TNUDG, 'Nudging time scale for sponge layers (days)', default=0.0) + call get_param(PF, mod, "ISOMIP_TNUDG", TNUDG, 'Nudging time scale for sponge layers (days)', default=0.0) - call get_param(PF, mod, "T_REF", t_ref, 'Reference temperature', default=10.0,& + call get_param(PF, mod, "ISOMIP_T_REF", t_ref, 'Reference temperature', default=10.0,& do_not_log=.true.) - call get_param(PF, mod, "S_REF", s_ref, 'Reference salinity', default=35.0,& + call get_param(PF, mod, "ISOMIP_S_REF", s_ref, 'Reference salinity', default=35.0,& do_not_log=.true.) - call get_param(PF, mod, "S_SUR_SPONGE", s_sur, 'Surface salinity in sponge layer.', default=s_ref) + call get_param(PF, mod, "ISOMIP_S_SUR_SPONGE", s_sur, 'Surface salinity in sponge layer.', default=s_ref) - call get_param(PF, mod, "S_BOT_SPONGE", s_bot, 'Bottom salinity in sponge layer.', default=s_ref) + call get_param(PF, mod, "ISOMIP_S_BOT_SPONGE", s_bot, 'Bottom salinity in sponge layer.', default=s_ref) - call get_param(PF, mod, "T_SUR_SPONGE", t_sur, 'Surface temperature in sponge layer.', default=t_ref) + call get_param(PF, mod, "ISOMIP_T_SUR_SPONGE", t_sur, 'Surface temperature in sponge layer.', default=t_ref) - call get_param(PF, mod, "T_BOT_SPONGE", t_bot, 'Bottom temperature in sponge layer.', default=t_ref) + call get_param(PF, mod, "ISOMIP_T_BOT_SPONGE", t_bot, 'Bottom temperature in sponge layer.', default=t_ref) T(:,:,:) = 0.0 ; S(:,:,:) = 0.0 ; Idamp(:,:) = 0.0; RHO(:,:,:) = 0.0 S_range = s_sur - s_bot @@ -378,11 +360,11 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) case ( REGRIDDING_LAYER, REGRIDDING_RHO ) ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT call calculate_density(t_sur,s_sur,0.0,rho_sur,tv%eqn_of_state) - write (*,*)'Surface density in sponge:', rho_sur + !write (*,*)'Surface density in sponge:', rho_sur call calculate_density(t_bot,s_bot,0.0,rho_bot,tv%eqn_of_state) - write (*,*)'Bottom density in sponge:', rho_bot + !write (*,*)'Bottom density in sponge:', rho_bot rho_range = rho_bot - rho_sur - write (*,*)'Density range in sponge:', rho_range + !write (*,*)'Density range in sponge:', rho_range ! Construct notional interface positions e0(1) = 0. @@ -487,7 +469,7 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) i=G%iec; j=G%jec do k = 1,nz call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state) - write(*,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) + !write(*,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) enddo case default @@ -501,9 +483,6 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) ! By default, momentum is advected vertically within the sponge, but ! ! momentum is typically not damped within the sponge. ! -! At this point, the ISOMIP configuration is done. The following are here as a -! template for other configurations. - ! The remaining calls to set_up_sponge_field can be in any order. ! if ( associated(tv%T) ) then call set_up_ALE_sponge_field(T,G,tv%T,CSp) From a7a1ede5800068c2a19ada1055a067eb081ab44c Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 27 May 2016 12:47:29 -0400 Subject: [PATCH 05/46] Removed ISOMIP_ from T_REF and S_REF --- src/user/ISOMIP_initialization.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 7c8264f74b..77b3ec2cd8 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -332,10 +332,10 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) call get_param(PF, mod, "ISOMIP_TNUDG", TNUDG, 'Nudging time scale for sponge layers (days)', default=0.0) - call get_param(PF, mod, "ISOMIP_T_REF", t_ref, 'Reference temperature', default=10.0,& + call get_param(PF, mod, "T_REF", t_ref, 'Reference temperature', default=10.0,& do_not_log=.true.) - call get_param(PF, mod, "ISOMIP_S_REF", s_ref, 'Reference salinity', default=35.0,& + call get_param(PF, mod, "S_REF", s_ref, 'Reference salinity', default=35.0,& do_not_log=.true.) call get_param(PF, mod, "ISOMIP_S_SUR_SPONGE", s_sur, 'Surface salinity in sponge layer.', default=s_ref) From 5708dbb46e2e3ef708deeaa508e85d30c2e4cba5 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 7 Jun 2016 16:57:06 -0400 Subject: [PATCH 06/46] Small changes in the initialization --- src/user/ISOMIP_initialization.F90 | 34 ++++++++++++++++++++++++------ 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 77b3ec2cd8..fe6455e1e8 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -113,10 +113,17 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) bx = 0.0; by = 0.0; xtil = 0.0 do j=js,je ; do i=is,ie + ! this is for 2D + !xtil = G%geoLonT(i,j)*1.0e3/xbar + !bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 + !by = (dc/(1.+exp(-2.*(1.*1.0e3- ly/2. - wc)/fc))) + & + ! (dc/(1.+exp(2.*(1.*1.0e3- ly/2. + wc)/fc))) + xtil = G%geoLonT(i,j)*1.0e3/xbar bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & - (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) + (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) + ! depth is positive D(i,j) = -max(bx+by,-bmax) @@ -264,7 +271,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, select case ( coordinateMode(verticalCoordinate) ) - case ( REGRIDDING_SIGMA, REGRIDDING_ZSTAR, REGRIDDING_RHO ) + case ( REGRIDDING_RHO, REGRIDDING_ZSTAR, REGRIDDING_SIGMA ) S_range = S_range / G%max_depth ! Convert S_range into dS/dz T_range = T_range / G%max_depth ! Convert T_range into dT/dz do j=js,je ; do i=is,ie @@ -276,11 +283,24 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, xi0 = xi0 + 0.5 * h(i,j,k) ! Depth at top of layer enddo enddo ; enddo -i=G%iec; j=G%jec -do k = 1,nz - call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,eqn_of_state) - !write(*,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) -enddo + +!i=G%iec; j=G%jec +!do k = 1,nz +! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,eqn_of_state) +! !write(*,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) +!enddo + +! case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA ) +! +! do j=js,je ; do i=is,ie +! xi0 = 0.0; +! do k = 1,nz +! xi1 = xi0 + h(i,j,k) / G%max_depth; +! S(i,j,k) = S_sur + 0.5 * S_range * (xi0 + xi1); +! T(i,j,k) = T_sur + 0.5 * T_range * (xi0 + xi1); +! xi0 = xi1; +! enddo +! enddo ; enddo case default call MOM_error(FATAL,"isomip_initialize: "// & From 4c12712132a79ab55fd2f9252cab605d659a04b6 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 9 Jun 2016 16:27:00 -0400 Subject: [PATCH 07/46] Added option for 2D verions of ISOMIP This option is very useful to develop/debug the model and can be activated with ISOMIP_2D = True. Note that the ice shelf .nc file also needs to be modified. --- src/user/ISOMIP_initialization.F90 | 49 ++++++++++++++++++------------ 1 file changed, 30 insertions(+), 19 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index fe6455e1e8..0d83cf4b06 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -86,7 +86,8 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) real :: wc ! half-width of the trough real :: ly ! domain width (across ice flow) real :: bx, by, xtil ! dummy vatiables - + logical :: is_2D ! If true, use 2D setup + ! G%ieg and G%jeg are the last indices in the global domain ! real, dimension (G%ieg) :: xtil ! eq. 3 ! real, dimension (G%ieg) :: bx ! eq. 2 @@ -106,30 +107,40 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) call log_version(param_file, mod, version, "") call get_param(param_file, mod, "MINIMUM_DEPTH", min_depth, & "The minimum depth of the ocean.", units="m", default=0.0) + call get_param(param_file, mod, "ISOMIP_2D",is_2D,'If true, use a 2D setup.', default=.false.) ! The following variables should be transformed into runtime parameters? bmax=720.0; b0=-150.0; b2=-728.8; b4=343.91; b6=-50.57 xbar=300.0E3; dc=500.0; fc=4.0E3; wc=24.0E3; ly=80.0E3 bx = 0.0; by = 0.0; xtil = 0.0 - do j=js,je ; do i=is,ie - ! this is for 2D - !xtil = G%geoLonT(i,j)*1.0e3/xbar - !bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 - !by = (dc/(1.+exp(-2.*(1.*1.0e3- ly/2. - wc)/fc))) + & - ! (dc/(1.+exp(2.*(1.*1.0e3- ly/2. + wc)/fc))) - - xtil = G%geoLonT(i,j)*1.0e3/xbar - bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 - by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & - (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) - -! depth is positive - D(i,j) = -max(bx+by,-bmax) - - if (D(i,j) > max_depth) D(i,j) = max_depth - if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth - enddo ; enddo + + if (is_2D) then + do j=js,je ; do i=is,ie + ! 2D setup + xtil = G%geoLonT(i,j)*1.0e3/xbar + bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 + by = (dc/(1.+exp(-2.*(40.*1.0e3- ly/2. - wc)/fc))) + & + (dc/(1.+exp(2.*(40.*1.0e3- ly/2. + wc)/fc))) + + D(i,j) = -max(bx+by,-bmax) + if (D(i,j) > max_depth) D(i,j) = max_depth + if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth + enddo ; enddo + + else + do j=js,je ; do i=is,ie + ! 3D + xtil = G%geoLonT(i,j)*1.0e3/xbar + bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 + by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & + (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) + + D(i,j) = -max(bx+by,-bmax) + if (D(i,j) > max_depth) D(i,j) = max_depth + if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth + enddo ; enddo + endif end subroutine ISOMIP_initialize_topography ! ----------------------------------------------------------------------------- From c1e088253199eadb5c3f74f497bd6dea5c91001d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 15 Jun 2016 10:31:02 -0400 Subject: [PATCH 08/46] 2D version of ISOMIP Just cleaned the code. --- src/user/ISOMIP_initialization.F90 | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 0d83cf4b06..ca0313984f 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -119,9 +119,13 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) do j=js,je ; do i=is,ie ! 2D setup xtil = G%geoLonT(i,j)*1.0e3/xbar + !xtil = 450*1.0e3/xbar bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 - by = (dc/(1.+exp(-2.*(40.*1.0e3- ly/2. - wc)/fc))) + & - (dc/(1.+exp(2.*(40.*1.0e3- ly/2. + wc)/fc))) + !by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & + ! (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) + + by = (dc/(1.+exp(-2.*(20.*1.0e3- ly/2. - wc)/fc))) + & + (dc/(1.+exp(2.*(20.*1.0e3- ly/2. + wc)/fc))) D(i,j) = -max(bx+by,-bmax) if (D(i,j) > max_depth) D(i,j) = max_depth @@ -131,7 +135,7 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) else do j=js,je ; do i=is,ie ! 3D - xtil = G%geoLonT(i,j)*1.0e3/xbar + !xtil = G%geoLonT(i,j)*1.0e3/xbar bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) From b2cd3d395f6fdec0273fed997de0396607bf48bc Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 16 Jun 2016 11:57:54 -0400 Subject: [PATCH 09/46] Sponge for layer mode in ISOMIP test case Changed the call in the ISOMIP_initialize_sponges to include the option to apply a sponge layer in layer modei (i.e., using sponge_CSp instead of ALE_sponge_CSp). --- src/initialization/MOM_state_initialization.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 08825eaff5..37bc9dafb6 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -409,7 +409,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & case ("DOME"); call DOME_initialize_sponges(G, GV, tv, PF, sponge_CSp) case ("DOME2D"); call DOME2d_initialize_sponges(G, GV, tv, PF, useALE, & sponge_CSp, ALE_sponge_CSp) - case ("ISOMIP"); call ISOMIP_initialize_sponges(G, GV, tv, PF, ALE_sponge_CSp) + case ("ISOMIP"); call ISOMIP_initialize_sponges(G, GV, tv, PF, useALE, & + sponge_CSp, ALE_sponge_CSp) case ("USER"); call user_initialize_sponges(G, use_temperature, tv, & PF, sponge_CSp, h) case ("phillips"); call Phillips_initialize_sponges(G, use_temperature, tv, & From 0c3a3a37e08f33fdfdfd24aa99cdb699c881ca71 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 16 Jun 2016 17:32:56 -0400 Subject: [PATCH 10/46] Sponge for layer mode in ISOMIP test case * Changed the call in the ISOMIP_initialize_sponges to include the option to appl y a sponge layer in layer modei (i.e., using sponge_CSp instead of ALE_sponge_CSp). * Changed sponge subroutine in ISOMIP_initialization to include the option for sponge when using layer mode. The model compiles, but sponges currently do not work for the ISOMIP in layer mode. --- src/user/ISOMIP_initialization.F90 | 249 ++++++++++++++++------------- 1 file changed, 137 insertions(+), 112 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index ca0313984f..3f158eb14b 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -25,6 +25,7 @@ module ISOMIP_initialization !* * !********+*********+*********+*********+*********+*********+*********+** use MOM_ALE_sponge, only : ALE_sponge_CS, set_up_ALE_sponge_field, initialize_ALE_sponge +use MOM_sponge, only : sponge_CS, set_up_sponge_field, initialize_sponge use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories @@ -89,10 +90,6 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) logical :: is_2D ! If true, use 2D setup ! G%ieg and G%jeg are the last indices in the global domain -! real, dimension (G%ieg) :: xtil ! eq. 3 -! real, dimension (G%ieg) :: bx ! eq. 2 -! real, dimension (G%ieg,G%jeg) :: dummy1 ! dummy array -! real, dimension (G%jeg) :: by ! eq. 4 ! This include declares and sets the variable "version". #include "version_variable.h" @@ -124,8 +121,9 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) !by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & ! (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) - by = (dc/(1.+exp(-2.*(20.*1.0e3- ly/2. - wc)/fc))) + & - (dc/(1.+exp(2.*(20.*1.0e3- ly/2. + wc)/fc))) + ! slice at y = 40 km + by = (dc/(1.+exp(-2.*(40.0*1.0e3- ly/2. - wc)/fc))) + & + (dc/(1.+exp(2.*(40.0*1.0e3- ly/2. + wc)/fc))) D(i,j) = -max(bx+by,-bmax) if (D(i,j) > max_depth) D(i,j) = max_depth @@ -134,8 +132,8 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) else do j=js,je ; do i=is,ie - ! 3D - !xtil = G%geoLonT(i,j)*1.0e3/xbar + ! 3D setup + xtil = G%geoLonT(i,j)*1.0e3/xbar bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) @@ -328,12 +326,14 @@ end subroutine ISOMIP_initialize_temperature_salinity !> Sets up the the inverse restoration time (Idamp), and ! ! the values towards which the interface heights and an arbitrary ! ! number of tracers should be restored within each sponge. -subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) +subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The container for vertical grid data type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available thermodynamic fields, including potential temperature and salinity or mixed layer density. Absent fields have NULL ptrs. type(param_file_type), intent(in) :: PF !< A structure indicating the open file to parse for model parameter values. - type(ALE_sponge_CS), pointer :: CSp !< A pointer that is set to point to the control structure for this module. + logical, intent(in) :: use_ALE !< If true, indicates model is in ALE mode + type(sponge_CS), pointer :: CSp !< Layer-mode sponge structure + type(ALE_sponge_CS), pointer :: ACSp !< ALE-mode sponge structure real :: T(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for temp real :: S(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for salt @@ -346,9 +346,11 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) real :: t_ref, s_ref ! reference T and S real :: rho_sur, rho_bot, rho_range, t_range, s_range - real :: e0(SZK_(G)) ! The resting interface heights, in m, usually ! + real :: e0(SZK_(G)+1) ! The resting interface heights, in m, usually ! ! negative because it is positive upward. ! real :: eta1D(SZK_(G)+1) ! Interface height relative to the sea surface ! + real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1) ! A temporary array for eta. + ! positive upward, in m. real :: min_depth, dummy1, z, delta_h real :: damp, rho_dummy, min_thickness, rho_tmp, xi0 @@ -389,107 +391,102 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) call get_param(PF, mod, "MINIMUM_DEPTH", min_depth, & "The minimum depth of the ocean.", units="m", default=0.0) - - select case ( coordinateMode(verticalCoordinate) ) - - case ( REGRIDDING_LAYER, REGRIDDING_RHO ) - ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT - call calculate_density(t_sur,s_sur,0.0,rho_sur,tv%eqn_of_state) - !write (*,*)'Surface density in sponge:', rho_sur - call calculate_density(t_bot,s_bot,0.0,rho_bot,tv%eqn_of_state) - !write (*,*)'Bottom density in sponge:', rho_bot - rho_range = rho_bot - rho_sur - !write (*,*)'Density range in sponge:', rho_range - - ! Construct notional interface positions - e0(1) = 0. - do K=2,nz - e0(k) = -G%max_depth * ( 0.5 * ( GV%Rlay(k-1) + GV%Rlay(k) ) - rho_sur ) / rho_range - e0(k) = min( 0., e0(k) ) ! Bound by surface - e0(k) = max( -G%max_depth, e0(k) ) ! Bound by possible deepest point in model - !write(*,*)'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k) - - enddo - e0(nz+1) = -G%max_depth - - ! Calculate thicknesses - do j=js,je ; do i=is,ie - eta1D(nz+1) = -1.0*G%bathyT(i,j) - do k=nz,1,-1 - eta1D(k) = e0(k) - if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_z)) then - eta1D(k) = eta1D(k+1) + GV%Angstrom_z - h(i,j,k) = GV%Angstrom_z - else - h(i,j,k) = eta1D(k) - eta1D(k+1) - endif - enddo - enddo ; enddo - - case ( REGRIDDING_ZSTAR ) ! Initial thicknesses for z coordinates - do j=js,je ; do i=is,ie - eta1D(nz+1) = -1.0*G%bathyT(i,j) - do k=nz,1,-1 - eta1D(k) = -G%max_depth * real(k-1) / real(nz) - if (eta1D(k) < (eta1D(k+1) + min_thickness)) then - eta1D(k) = eta1D(k+1) + min_thickness - h(i,j,k) = min_thickness - else - h(i,j,k) = eta1D(k) - eta1D(k+1) - endif - enddo - enddo ; enddo - - case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates - do j=js,je ; do i=is,ie - delta_h = G%bathyT(i,j) / dfloat(nz) - h(i,j,:) = delta_h - end do ; end do + if (associated(CSp)) call MOM_error(FATAL, & + "ISOMIP_initialize_sponges called with an associated control structure.") + if (associated(ACSp)) call MOM_error(FATAL, & + "ISOMIP_initialize_sponges called with an associated ALE-sponge control structure.") - case default - call MOM_error(FATAL,"ISOMIP_initialize_sponges: "// & - "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE") - - end select - -! Here the inverse damping time, in s-1, is set. Set Idamp to 0 ! -! wherever there is no sponge, and the subroutines that are called ! -! will automatically set up the sponges only where Idamp is positive! -! and mask2dT is 1. + ! Here the inverse damping time, in s-1, is set. Set Idamp to 0 ! + ! wherever there is no sponge, and the subroutines that are called ! + ! will automatically set up the sponges only where Idamp is positive! + ! and mask2dT is 1. do i=is,ie; do j=js,je if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then -! 1 / day + ! 1 / day dummy1=(G%geoLonT(i,j)-790.0)/(800.0-790.0) damp = 1.0/TNUDG * max(0.0,dummy1) else ; damp=0.0 endif -! convert to 1 / seconds + ! convert to 1 / seconds if (G%bathyT(i,j) > min_depth) then Idamp(i,j) = damp/86400.0 else ; Idamp(i,j) = 0.0 ; endif - + enddo ; enddo - if (associated(CSp)) then - call MOM_error(WARNING, "ISOMIP_initialize_sponges called with an associated "// & - "control structure.") - return - endif + ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT + call calculate_density(t_sur,s_sur,0.0,rho_sur,tv%eqn_of_state) + !write (*,*)'Surface density in sponge:', rho_sur + call calculate_density(t_bot,s_bot,0.0,rho_bot,tv%eqn_of_state) + !write (*,*)'Bottom density in sponge:', rho_bot + rho_range = rho_bot - rho_sur + !write (*,*)'Density range in sponge:', rho_range + + + if (use_ALE) then + select case ( coordinateMode(verticalCoordinate) ) + + case ( REGRIDDING_RHO ) + ! Construct notional interface positions + e0(1) = 0. + do K=2,nz + e0(k) = -G%max_depth * ( 0.5 * ( GV%Rlay(k-1) + GV%Rlay(k) ) - rho_sur ) / rho_range + e0(k) = min( 0., e0(k) ) ! Bound by surface + e0(k) = max( -G%max_depth, e0(k) ) ! Bound by possible deepest point in model + !write(*,*)'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k) + + enddo + e0(nz+1) = -G%max_depth + + ! Calculate thicknesses + do j=js,je ; do i=is,ie + eta1D(nz+1) = -1.0*G%bathyT(i,j) + do k=nz,1,-1 + eta1D(k) = e0(k) + if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_z)) then + eta1D(k) = eta1D(k+1) + GV%Angstrom_z + h(i,j,k) = GV%Angstrom_z + else + h(i,j,k) = eta1D(k) - eta1D(k+1) + endif + enddo + enddo ; enddo + + case ( REGRIDDING_ZSTAR ) ! Initial thicknesses for z coordinates + do j=js,je ; do i=is,ie + eta1D(nz+1) = -1.0*G%bathyT(i,j) + do k=nz,1,-1 + eta1D(k) = -G%max_depth * real(k-1) / real(nz) + if (eta1D(k) < (eta1D(k+1) + min_thickness)) then + eta1D(k) = eta1D(k+1) + min_thickness + h(i,j,k) = min_thickness + else + h(i,j,k) = eta1D(k) - eta1D(k+1) + endif + enddo + enddo ; enddo -! This call sets up the damping rates and interface heights. -! This sets the inverse damping timescale fields in the sponges. ! + case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates + do j=js,je ; do i=is,ie + delta_h = G%bathyT(i,j) / dfloat(nz) + h(i,j,:) = delta_h + end do ; end do - call initialize_ALE_sponge(Idamp,h, nz, G, PF, CSp) + case default + call MOM_error(FATAL,"ISOMIP_initialize_sponges: "// & + "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE") -! setup temp and salt at the sponge zone - select case ( coordinateMode(verticalCoordinate) ) + end select + + ! This call sets up the damping rates and interface heights. + ! This sets the inverse damping timescale fields in the sponges. + call initialize_ALE_sponge(Idamp,h, nz, G, PF, ACSp) - case ( REGRIDDING_SIGMA, REGRIDDING_ZSTAR, REGRIDDING_RHO ) S_range = S_range / G%max_depth ! Convert S_range into dS/dz T_range = T_range / G%max_depth ! Convert T_range into dT/dz do j=js,je ; do i=is,ie @@ -501,33 +498,61 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) xi0 = xi0 + 0.5 * h(i,j,k) ! Depth at top of layer enddo enddo ; enddo - i=G%iec; j=G%jec - do k = 1,nz - call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state) - !write(*,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) - enddo - - case default - call MOM_error(FATAL,"isomip_initialize_sponge: "// & - "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE") - - end select + ! for debugging + !i=G%iec; j=G%jec + !do k = 1,nz + ! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state) + !write(*,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) + !enddo + + ! Now register all of the fields which are damped in the sponge. ! + ! By default, momentum is advected vertically within the sponge, but ! + ! momentum is typically not damped within the sponge. ! + + ! The remaining calls to set_up_sponge_field can be in any order. ! + if ( associated(tv%T) ) then + call set_up_ALE_sponge_field(T,G,tv%T,ACSp) + endif + if ( associated(tv%S) ) then + call set_up_ALE_sponge_field(S,G,tv%S,ACSp) + endif -! Now register all of the fields which are damped in the sponge. ! -! By default, momentum is advected vertically within the sponge, but ! -! momentum is typically not damped within the sponge. ! -! The remaining calls to set_up_sponge_field can be in any order. ! - if ( associated(tv%T) ) then - call set_up_ALE_sponge_field(T,G,tv%T,CSp) - endif + else ! layer mode + + ! Construct notional interface positions + e0(1) = 0. + do K=2,nz + e0(k) = -G%max_depth * ( 0.5 * ( GV%Rlay(k-1) + GV%Rlay(k) ) - rho_sur ) / rho_range + e0(k) = min( 0., e0(k) ) ! Bound by surface + e0(k) = max( -G%max_depth, e0(k) ) ! Bound by possible deepest point in model + !write(*,*)'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k) + enddo + e0(nz+1) = -G%max_depth + + ! Calculate thicknesses + do j=js,je ; do i=is,ie + eta1D(nz+1) = -1.0*G%bathyT(i,j) + do k=nz,1,-1 + eta1D(k) = e0(k) + if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_z)) then + eta1D(k) = eta1D(k+1) + GV%Angstrom_z + h(i,j,k) = GV%Angstrom_z + else + h(i,j,k) = eta1D(k) - eta1D(k+1) + endif + enddo + + eta(i,j,nz+1) = -G%bathyT(i,j) + do K=nz,1,-1 + eta(i,j,K) = eta(i,j,K+1) + h(i,j,k) + enddo + enddo ; enddo + call initialize_sponge(Idamp, eta, G, PF, CSp) - if ( associated(tv%S) ) then - call set_up_ALE_sponge_field(S,G,tv%S,CSp) endif - end subroutine ISOMIP_initialize_sponges ! ----------------------------------------------------------------------------- From 5257667e1147f13a2bb8d0f201f76dc3fac32812 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 17 Jun 2016 11:20:55 -0400 Subject: [PATCH 11/46] Change G to HI in ISOMIP tracer * Changed ISOMIP_tracer.F90 and ISOMIP call in MOM_tracer_flow_control.F90 following this (this commit)[https://github.com/NOAA-GFDL/MOM6/commit/401b370921461bcb5c5f023284678d279ee3236b] * Fix conflit in ISOMIP_initialization.F90 Model compiles sucessfully! --- src/tracer/ISOMIP_tracer.F90 | 12 ++--- src/tracer/MOM_tracer_flow_control.F90 | 8 +-- src/user/ISOMIP_initialization.F90 | 71 +++++++++++++------------- 3 files changed, 46 insertions(+), 45 deletions(-) diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index 152980b04b..072f5b4d87 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -21,6 +21,7 @@ module ISOMIP_tracer use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_forcing_type, only : forcing +use MOM_hor_index, only : hor_index_type use MOM_grid, only : ocean_grid_type use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc use MOM_restart, only : register_restart_field, MOM_restart_CS @@ -87,12 +88,12 @@ module ISOMIP_tracer !> This subroutine is used to register tracer fields -function register_ISOMIP_tracer(G, param_file, CS, diag, tr_Reg, & +function register_ISOMIP_tracer(HI, GV, param_file, CS, tr_Reg, & restart_CS) - type(ocean_grid_type), intent(in) :: G ! NULL() logical :: register_ISOMIP_tracer integer :: isd, ied, jsd, jed, nz, m - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = G%ke + isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke if (associated(CS)) then call MOM_error(WARNING, "ISOMIP_register_tracer called with an "// & @@ -113,7 +114,6 @@ function register_ISOMIP_tracer(G, param_file, CS, diag, tr_Reg, & endif allocate(CS) - CS%diag => diag ! Read all relevant parameters and write them to the model log. call log_version(param_file, mod, version, "") call get_param(param_file, mod, "ISOMIP_TRACER_IC_FILE", CS%tracer_IC_file, & @@ -149,7 +149,7 @@ function register_ISOMIP_tracer(G, param_file, CS, diag, tr_Reg, & ! Register the tracer for the restart file. call register_restart_field(tr_ptr, CS%tr_desc(m), .true., restart_CS) ! Register the tracer for horizontal advection & diffusion. - call register_tracer(tr_ptr, CS%tr_desc(m), param_file, G, tr_Reg, & + call register_tracer(tr_ptr, CS%tr_desc(m), param_file, HI, GV, tr_Reg, & tr_desc_ptr=CS%tr_desc(m)) ! Set coupled_tracers to be true (hard-coded above) to provide the surface diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index d6fb9c4e7c..3b1591d031 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -184,11 +184,11 @@ subroutine call_tracer_register(HI, GV, param_file, CS, tr_Reg, restart_CS) USER_register_tracer_example(HI, GV, param_file, CS%USER_tracer_example_CSp, & tr_Reg, restart_CS) if (CS%use_DOME_tracer) CS%use_DOME_tracer = & - register_DOME_tracer(G, param_file, CS%DOME_tracer_CSp, & - diag, tr_Reg, restart_CS) + register_DOME_tracer(HI, GV, param_file, CS%DOME_tracer_CSp, & + tr_Reg, restart_CS) if (CS%use_ISOMIP_tracer) CS%use_ISOMIP_tracer = & - register_ISOMIP_tracer(G, param_file, CS%ISOMIP_tracer_CSp, & - diag, tr_Reg, restart_CS) + register_ISOMIP_tracer(HI, GV, param_file, CS%ISOMIP_tracer_CSp, & + tr_Reg, restart_CS) if (CS%use_ideal_age) CS%use_ideal_age = & register_ideal_age_tracer(HI, GV, param_file, CS%ideal_age_tracer_CSp, & tr_Reg, restart_CS) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index acda496889..bfa4242ca8 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -432,8 +432,8 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) rho_range = rho_bot - rho_sur !write (*,*)'Density range in sponge:', rho_range - - if (use_ALE) then + if (use_ALE) then + select case ( coordinateMode(verticalCoordinate) ) case ( REGRIDDING_RHO ) @@ -482,43 +482,44 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) h(i,j,:) = delta_h end do ; end do - case default + case default call MOM_error(FATAL,"ISOMIP_initialize_sponges: "// & "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE") - ! This call sets up the damping rates and interface heights. - ! This sets the inverse damping timescale fields in the sponges. - call initialize_ALE_sponge(Idamp,h, nz, G, PF, ACSp) + end select + ! This call sets up the damping rates and interface heights. + ! This sets the inverse damping timescale fields in the sponges. + call initialize_ALE_sponge(Idamp,h, nz, G, PF, ACSp) - S_range = S_range / G%max_depth ! Convert S_range into dS/dz - T_range = T_range / G%max_depth ! Convert T_range into dT/dz - do j=js,je ; do i=is,ie - xi0 = -G%bathyT(i,j); - do k = nz,1,-1 - xi0 = xi0 + 0.5 * h(i,j,k) ! Depth in middle of layer - S(i,j,k) = S_sur + S_range * xi0 - T(i,j,k) = T_sur + T_range * xi0 - xi0 = xi0 + 0.5 * h(i,j,k) ! Depth at top of layer - enddo - enddo ; enddo - ! for debugging - !i=G%iec; j=G%jec - !do k = 1,nz - ! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state) - !write(*,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) - !enddo - - ! Now register all of the fields which are damped in the sponge. ! - ! By default, momentum is advected vertically within the sponge, but ! - ! momentum is typically not damped within the sponge. ! - - ! The remaining calls to set_up_sponge_field can be in any order. ! - if ( associated(tv%T) ) then - call set_up_ALE_sponge_field(T,G,tv%T,ACSp) - endif - if ( associated(tv%S) ) then - call set_up_ALE_sponge_field(S,G,tv%S,ACSp) - endif + S_range = S_range / G%max_depth ! Convert S_range into dS/dz + T_range = T_range / G%max_depth ! Convert T_range into dT/dz + do j=js,je ; do i=is,ie + xi0 = -G%bathyT(i,j); + do k = nz,1,-1 + xi0 = xi0 + 0.5 * h(i,j,k) ! Depth in middle of layer + S(i,j,k) = S_sur + S_range * xi0 + T(i,j,k) = T_sur + T_range * xi0 + xi0 = xi0 + 0.5 * h(i,j,k) ! Depth at top of layer + enddo + enddo ; enddo + ! for debugging + !i=G%iec; j=G%jec + !do k = 1,nz + ! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state) + !write(*,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) + !enddo + + ! Now register all of the fields which are damped in the sponge. ! + ! By default, momentum is advected vertically within the sponge, but ! + ! momentum is typically not damped within the sponge. ! + + ! The remaining calls to set_up_sponge_field can be in any order. ! + if ( associated(tv%T) ) then + call set_up_ALE_sponge_field(T,G,tv%T,ACSp) + endif + if ( associated(tv%S) ) then + call set_up_ALE_sponge_field(S,G,tv%S,ACSp) + endif else ! layer mode From 2c9b75eb35fba2fc3b134d207699075116805033 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 21 Jun 2016 17:16:37 -0400 Subject: [PATCH 12/46] Allocated additional fluxes in MOM_ice_shelf To use ENERGETICS_SFC_PBL = true, additional fluxes had to be allocated in MOM_ice_shelf.F90. This was done using the function allocate_forcing_type. --- src/ice_shelf/MOM_ice_shelf.F90 | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index ef2b752cd9..57535beb44 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -862,7 +862,6 @@ subroutine add_shelf_flux(G, CS, state, fluxes) tauy2 = (asv1 * state%tauy_shelf(i,j-1)**2 + & asv2 * state%tauy_shelf(i,j)**2 ) / (asv1 + asv2) fluxes%ustar(i,j) = MAX(CS%ustar_bg, sqrt(Irho0 * sqrt(taux2 + tauy2))) - if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0 if (associated(fluxes%lw)) fluxes%lw(i,j) = 0.0 if (associated(fluxes%latent)) fluxes%latent(i,j) = 0.0 @@ -1380,16 +1379,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti ! Allocate the arrays for passing ice-shelf data through the forcing type. if (.not. solo_ice_sheet) then if (is_root_pe()) print *,"initialize_ice_shelf: allocating fluxes" - allocate( fluxes%frac_shelf_h(isd:ied,jsd:jed) ) ; fluxes%frac_shelf_h(:,:) = 0.0 - allocate( fluxes%frac_shelf_u(Isdq:Iedq,jsd:jed) ) ; fluxes%frac_shelf_u(:,:) = 0.0 - allocate( fluxes%frac_shelf_v(isd:ied,Jsdq:Jedq) ) ; fluxes%frac_shelf_v(:,:) = 0.0 - allocate( fluxes%ustar_shelf(isd:ied,jsd:jed) ) ; fluxes%ustar_shelf(:,:) = 0.0 - allocate( fluxes%iceshelf_melt(isd:ied,jsd:jed) ) ; fluxes%iceshelf_melt(:,:) = 0.0 - if (.not.associated(fluxes%p_surf)) then - allocate( fluxes%p_surf(isd:ied,jsd:jed) ) ; fluxes%p_surf(:,:) = 0.0 - endif - allocate( fluxes%rigidity_ice_u(Isdq:Iedq,jsd:jed) ) ; fluxes%rigidity_ice_u(:,:) = 0.0 - allocate( fluxes%rigidity_ice_v(isd:ied,Jsdq:Jedq) ) ; fluxes%rigidity_ice_v(:,:) = 0.0 + !GM, the following is necessary to make ENERGETICS_SFC_PBL work + call allocate_forcing_type(G, fluxes, ustar=.true., shelf=.true., & + press=.true., water=.true., heat=.true.) else if (is_root_pe()) print *,"allocating fluxes in solo mode" call allocate_forcing_type(G, fluxes, ustar=.true., shelf=.true., press=.true.) From 9e1453069311ad49a82aa2b684c6c5b179bcc1e4 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 21 Jun 2016 19:32:49 -0400 Subject: [PATCH 13/46] Added diag to the call and fixed missing pointer assignment diag was not been passed to initialize_ISOMIP_tracer. It needs to be passed so that the structure CS%diag can be used during the register_diag_field calls. --- src/tracer/ISOMIP_tracer.F90 | 9 ++++++--- src/tracer/MOM_tracer_flow_control.F90 | 2 +- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index 072f5b4d87..e63037286e 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -166,13 +166,15 @@ end function register_ISOMIP_tracer !> Initializes the NTR tracer fields in tr(:,:,:,:) ! and it sets up the tracer output. -subroutine initialize_ISOMIP_tracer(restart, day, G, GV, h, OBC, CS, ALE_sponge_CSp, & - diag_to_Z_CSp) +subroutine initialize_ISOMIP_tracer(restart, day, G, GV, h, diag, OBC, CS, & + ALE_sponge_CSp, diag_to_Z_CSp) + type(ocean_grid_type), intent(in) :: G !< Grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. logical, intent(in) :: restart !< .true. if the fields have already been read from a restart file. type(time_type), target, intent(in) :: day !< Time of the start of the run. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(diag_ctrl), target, intent(in) :: diag type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, where, and what open boundary conditions are used. This is not being used for now. type(ISOMIP_tracer_CS), pointer :: CS !< The control structure returned by a previous call to ISOMIP_register_tracer. type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< A pointer to the control structure for the sponges, if they are in use. Otherwise this may be unassociated. @@ -206,6 +208,7 @@ subroutine initialize_ISOMIP_tracer(restart, day, G, GV, h, OBC, CS, ALE_sponge_ h_neglect = GV%H_subroundoff CS%Time => day + CS%diag => diag if (.not.restart) then if (len_trim(CS%tracer_IC_file) >= 1) then @@ -340,7 +343,7 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G if (fluxes%iceshelf_melt(i,j) > 0.0) then ! CS%tr(i,j,1,m) = (dt*fluxes%iceshelf_melt(i,j)/(365.0 * 86400.0) & ! + h_old(i,j,1)*CS%tr(i,j,1,m))/h_new(i,j,1) - CS%tr(i,j,1,m) = 1.0 + CS%tr(i,j,1:3,m) = 1.0 endif enddo ; enddo enddo diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index 3b1591d031..f8299cf0a5 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -256,7 +256,7 @@ subroutine tracer_flow_control_init(restart, day, G, GV, h, param_file, diag, OB call initialize_DOME_tracer(restart, day, G, GV, h, diag, OBC, CS%DOME_tracer_CSp, & sponge_CSp, diag_to_Z_CSp) if (CS%use_ISOMIP_tracer) & - call initialize_ISOMIP_tracer(restart, day, G, GV, h, OBC, CS%ISOMIP_tracer_CSp, & + call initialize_ISOMIP_tracer(restart, day, G, GV, h, diag, OBC, CS%ISOMIP_tracer_CSp, & ALE_sponge_CSp, diag_to_Z_CSp) if (CS%use_ideal_age) & call initialize_ideal_age_tracer(restart, day, G, GV, h, diag, OBC, CS%ideal_age_tracer_CSp, & From e09aa158abad37a9e51f853282de6e6d98d1ff53 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 24 Jun 2016 16:53:53 -0400 Subject: [PATCH 14/46] Added sponge layers when using layer mode Now the ISOMIP initialization allows the implementation of sponge layer when using layer mode. This is done by reading a NETCDF file, for example the IC file saved by the model, that specifies eta, T and S. The T and S values most be consistent with the target density of the layers. --- src/user/ISOMIP_initialization.F90 | 77 +++++++++++++++++------------- 1 file changed, 45 insertions(+), 32 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index bfa4242ca8..1d943fe1be 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -359,7 +359,8 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) ! positive upward, in m. real :: min_depth, dummy1, z, delta_h real :: damp, rho_dummy, min_thickness, rho_tmp, xi0 - character(len=40) :: verticalCoordinate + character(len=40) :: verticalCoordinate, filename, state_file, inputdir + character(len=40) :: temp_var, salt_var, eta_var character(len=40) :: mod = "ISOMIP_initialize_sponges" ! This subroutine's name. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz @@ -521,37 +522,49 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) call set_up_ALE_sponge_field(S,G,tv%S,ACSp) endif - else ! layer mode - - ! Construct notional interface positions - e0(1) = 0. - do K=2,nz - e0(k) = -G%max_depth * ( 0.5 * ( GV%Rlay(k-1) + GV%Rlay(k) ) - rho_sur ) / rho_range - e0(k) = min( 0., e0(k) ) ! Bound by surface - e0(k) = max( -G%max_depth, e0(k) ) ! Bound by possible deepest point in model - !write(*,*)'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k) - enddo - e0(nz+1) = -G%max_depth - - ! Calculate thicknesses - do j=js,je ; do i=is,ie - eta1D(nz+1) = -1.0*G%bathyT(i,j) - do k=nz,1,-1 - eta1D(k) = e0(k) - if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_z)) then - eta1D(k) = eta1D(k+1) + GV%Angstrom_z - h(i,j,k) = GV%Angstrom_z - else - h(i,j,k) = eta1D(k) - eta1D(k+1) - endif - enddo - - eta(i,j,nz+1) = -G%bathyT(i,j) - do K=nz,1,-1 - eta(i,j,K) = eta(i,j,K+1) + h(i,j,k) - enddo - enddo ; enddo - call initialize_sponge(Idamp, eta, G, PF, CSp) + else ! layer mode + ! 1) Read eta, salt and temp from IC file + call get_param(PF, mod, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + call get_param(PF, mod, "ISOMIP_SPONGE_FILE", state_file, & + "The name of the file with the state to damp toward.", & + fail_if_missing=.true.) + write(*,*)'state_file', state_file + + call get_param(PF, mod, "SPONGE_PTEMP_VAR", temp_var, & + "The name of the potential temperature variable in \n"//& + "SPONGE_STATE_FILE.", default="Temp") + call get_param(PF, mod, "SPONGE_SALT_VAR", salt_var, & + "The name of the salinity variable in \n"//& + "SPONGE_STATE_FILE.", default="Salt") + call get_param(PF, mod, "SPONGE_ETA_VAR", eta_var, & + "The name of the interface height variable in \n"//& + "SPONGE_STATE_FILE.", default="eta") + + filename = trim(inputdir)//trim(state_file) + + if (.not.file_exists(filename, G%Domain)) & + call MOM_error(FATAL, " ISOMIP_initialize_sponges: Unable to open "//trim(filename)) + + !call read_data("",eta_var,eta(:,:,:), domain=G%Domain%mpp_domain) + call read_data(filename,eta_var,eta(:,:,:), domain=G%Domain%mpp_domain) + call read_data(filename,temp_var,T(:,:,:), domain=G%Domain%mpp_domain) + call read_data(filename,salt_var,S(:,:,:), domain=G%Domain%mpp_domain) + + ! for debugging + !i=G%iec; j=G%jec + !do k = 1,nz + ! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state) + ! write(*,*) 'Sponge - k,eta,T,S,rho,Rlay',k,eta(i,j,k),T(i,j,k),& + ! S(i,j,k),rho_tmp,GV%Rlay(k) + !enddo + + ! Set the inverse damping rates so that the model will know where to + ! apply the sponges, along with the interface heights. + call initialize_sponge(Idamp, eta, G, PF, CSp) + ! Apply sponge in tracer fields + call set_up_sponge_field(T, tv%T, G, nz, CSp) + call set_up_sponge_field(S, tv%S, G, nz, CSp) endif From 2715e54c480400eaae48bdfdafc6083cc95255f1 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 29 Jun 2016 09:27:57 -0400 Subject: [PATCH 15/46] Added sponge layer for layer mode In the ISOMIP layer mode, T/S is now restored within a sponge layer. This is done by first reading a netcdf file with T and S, then inverting to T or S (depending on FIT_SALINITY) to match the layer density. --- src/user/ISOMIP_initialization.F90 | 133 ++++++++++++++++++++++------- 1 file changed, 100 insertions(+), 33 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 1d943fe1be..ea6e9ea9c6 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -125,8 +125,17 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) else do j=js,je ; do i=is,ie - ! 3D setup + ! 3D setup + ! #### TEST ####### + !if (G%geoLonT(i,j)<500.) then + ! xtil = 500.*1.0e3/xbar + !else + ! xtil = G%geoLonT(i,j)*1.0e3/xbar + !endif + ! ##### TEST ##### + xtil = G%geoLonT(i,j)*1.0e3/xbar + bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) @@ -252,37 +261,41 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, type(EOS_type), pointer :: eqn_of_state !< Equation of state structure ! Local variables - integer :: i, j, k, is, ie, js, je, nz + integer :: i, j, k, is, ie, js, je, nz, itt real :: x, ds, dt, rho_sur, rho_bot real :: xi0, xi1, dxi, r, S_sur, T_sur, S_bot, T_bot, S_range, T_range real :: z ! vertical position in z space character(len=40) :: verticalCoordinate, density_profile real :: rho_tmp - + logical :: fit_salin ! If true, accept the prescribed temperature and fit the salinity. + real :: T0(SZK_(G)), S0(SZK_(G)) + real :: drho_dT(SZK_(G)) ! Derivative of density with temperature in kg m-3 K-1. ! + real :: drho_dS(SZK_(G)) ! Derivative of density with salinity in kg m-3 PSU-1. ! + real :: rho_guess(SZK_(G)) ! Potential density at T0 & S0 in kg m-3. + real :: pres(SZK_(G)) ! An array of the reference pressure in Pa. (zero here) + real :: drho_dT1, drho_dS1 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + pres(:) = 0.0 call get_param(param_file,mod,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & default=DEFAULT_COORDINATE_MODE) - call get_param(param_file, mod, "ISOMIP_T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9,& - do_not_log=.true.) - call get_param(param_file, mod, "ISOMIP_S_SUR", s_sur, 'Salinity at the surface (interface)', default=33.8,& - do_not_log=.true.) - call get_param(param_file, mod, "ISOMIP_T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9,& - do_not_log=.true.) - call get_param(param_file, mod, "ISOMIP_S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55,& - do_not_log=.true.) + call get_param(param_file, mod, "ISOMIP_T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9) + call get_param(param_file, mod, "ISOMIP_S_SUR", s_sur, 'Salinity at the surface (interface)', default=33.8) + call get_param(param_file, mod, "ISOMIP_T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9) + call get_param(param_file, mod, "ISOMIP_S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55) call calculate_density(t_sur,s_sur,0.0,rho_sur,eqn_of_state) !write (*,*)'Density in the surface layer:', rho_sur call calculate_density(t_bot,s_bot,0.0,rho_bot,eqn_of_state) !write (*,*)'Density in the bottom layer::', rho_bot - S_range = s_sur - s_bot - T_range = t_sur - t_bot - select case ( coordinateMode(verticalCoordinate) ) case ( REGRIDDING_RHO, REGRIDDING_ZSTAR, REGRIDDING_SIGMA ) + S_range = s_sur - s_bot + T_range = t_sur - t_bot + !write(*,*)'S_range,T_range',S_range,T_range + S_range = S_range / G%max_depth ! Convert S_range into dS/dz T_range = T_range / G%max_depth ! Convert T_range into dT/dz do j=js,je ; do i=is,ie @@ -295,29 +308,84 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, enddo enddo ; enddo -!i=G%iec; j=G%jec -!do k = 1,nz -! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,eqn_of_state) -! !write(*,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) -!enddo - -! case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA ) -! -! do j=js,je ; do i=is,ie -! xi0 = 0.0; -! do k = 1,nz -! xi1 = xi0 + h(i,j,k) / G%max_depth; -! S(i,j,k) = S_sur + 0.5 * S_range * (xi0 + xi1); -! T(i,j,k) = T_sur + 0.5 * T_range * (xi0 + xi1); -! xi0 = xi1; -! enddo -! enddo ; enddo - + case ( REGRIDDING_LAYER ) + call get_param(param_file, mod, "FIT_SALINITY", fit_salin, & + "If true, accept the prescribed temperature and fit the \n"//& + "salinity; otherwise take salinity and fit temperature.", & + default=.false.) + call get_param(param_file, mod, "DRHO_DS", drho_dS1, & + "Partial derivative of density with salinity.", & + units="kg m-3 PSU-1", fail_if_missing=.true.) + call get_param(param_file, mod, "DRHO_DT", drho_dT1, & + "Partial derivative of density with temperature.", & + units="kg m-3 K-1", fail_if_missing=.true.) + + !write(*,*)'read drho_dS, drho_dT', drho_dS1, drho_dT1 + + S_range = s_bot - s_sur + T_range = t_bot - t_sur + !write(*,*)'S_range,T_range',S_range,T_range + S_range = S_range / G%max_depth ! Convert S_range into dS/dz + T_range = T_range / G%max_depth ! Convert T_range into dT/dz + + i=G%iec; j=G%jec; xi0 = 0.0; + do k = 1,nz + xi1 = xi0 + 0.5 * h(i,j,k); + S0(k) = S_sur + S_range * xi1; + T0(k) = T_sur + T_range * xi1; + xi0 = xi0 + h(i,j,k); + !write(*,*)'S,T,xi0,xi1,k',S0(k),T0(k),xi0,xi1,k + enddo + + call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,1,eqn_of_state) + !write(*,*)'computed drho_dS, drho_dT', drho_dS(1), drho_dT(1) + call calculate_density(T0(1),S0(1),0.,rho_guess(1),eqn_of_state) + + if (fit_salin) then + ! A first guess of the layers' salinity. + do k=nz,1,-1 + S0(k) = max(0.0, S0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dS1) + enddo + ! Refine the guesses for each layer. + do itt=1,6 + call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state) + call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state) + do k=1,nz + S0(k) = max(0.0, S0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dS1) + enddo + enddo + + else + ! A first guess of the layers' temperatures. + do k=nz,1,-1 + T0(k) = T0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dT1 + enddo + + do itt=1,6 + call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state) + call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state) + do k=1,nz + T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) + enddo + enddo + endif + + do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + T(i,j,k) = T0(k) ; S(i,j,k) = S0(k) + enddo ; enddo ; enddo + case default call MOM_error(FATAL,"isomip_initialize: "// & "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE") end select + + ! for debugging + !i=G%iec; j=G%jec + !do k = 1,nz + ! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,eqn_of_state) + ! write(*,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) + !enddo end subroutine ISOMIP_initialize_temperature_salinity @@ -529,7 +597,6 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) call get_param(PF, mod, "ISOMIP_SPONGE_FILE", state_file, & "The name of the file with the state to damp toward.", & fail_if_missing=.true.) - write(*,*)'state_file', state_file call get_param(PF, mod, "SPONGE_PTEMP_VAR", temp_var, & "The name of the potential temperature variable in \n"//& From e9783ea4ba120e10ad77f3f850620bbe34fc4b7d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 30 Jun 2016 08:24:54 -0400 Subject: [PATCH 16/46] Allocated fluxes%buoy in MOM_ice_shelf to avoid [this](https://github.com/NOAA-GFDL/MOM6/blob/dev/master/src/parameterizations/vertical/MOM_entrain_diffusive.F90#L277) fatal error when using pure layer mode without defining BULKMIXEDLAYER. I also modified MOM_ALE_sponge.F90 so it reads REMAPPING_SCHEME from the input file to keep remapping consistent with the interior. --- src/ice_shelf/MOM_ice_shelf.F90 | 3 +++ src/parameterizations/vertical/MOM_ALE_sponge.F90 | 9 +++++++-- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 57535beb44..d57c847036 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1382,6 +1382,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti !GM, the following is necessary to make ENERGETICS_SFC_PBL work call allocate_forcing_type(G, fluxes, ustar=.true., shelf=.true., & press=.true., water=.true., heat=.true.) + ! This is to make it work with layer mode without BULKMIXEDLAYER + allocate( fluxes%buoy(isd:ied,jsd:jed) ) + fluxes%buoy(:,:) = 0.0 else if (is_root_pe()) print *,"allocating fluxes in solo mode" call allocate_forcing_type(G, fluxes, ustar=.true., shelf=.true., press=.true.) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index ce5f0778be..c56d4563a3 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -85,7 +85,7 @@ subroutine initialize_ALE_sponge(Iresttime, data_h, nz_data, G, param_file, CS) character(len=40) :: mod = "MOM_sponge" ! This module's name. logical :: use_sponge integer :: i, j, k, col, total_sponge_cols - + character(len=10) :: remapScheme if (associated(CS)) then call MOM_error(WARNING, "initialize_sponge called with an associated "// & "control structure.") @@ -103,6 +103,11 @@ subroutine initialize_ALE_sponge(Iresttime, data_h, nz_data, G, param_file, CS) allocate(CS) + call get_param(param_file, mod, "REMAPPING_SCHEME", remapScheme, & + "This sets the reconstruction scheme used \n"//& + " for vertical remapping for all variables.", & + default="PLM", do_not_log=.true.) + CS%nz = G%ke CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed @@ -144,7 +149,7 @@ subroutine initialize_ALE_sponge(Iresttime, data_h, nz_data, G, param_file, CS) call sum_across_PEs(total_sponge_cols) ! Call the constructor for remapping control structure - call initialize_remapping(CS%remap_cs, 'PPM_H4', boundary_extrapolation=.false.) + call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=.false.) call log_param(param_file, mod, "!Total sponge columns", total_sponge_cols, & "The total number of columns where sponges are applied.") From 9e9cad2829f3a15ba4ab0551fd881707b7edfdae Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 30 Jun 2016 11:43:22 -0400 Subject: [PATCH 17/46] Some changes that could not be tested because the model is broken. --- src/ice_shelf/MOM_ice_shelf.F90 | 8 +++++--- src/user/ISOMIP_initialization.F90 | 1 + 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 9af35dd490..c488800ce7 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1389,11 +1389,13 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti if (.not. solo_ice_sheet) then if (is_root_pe()) print *,"initialize_ice_shelf: allocating fluxes" !GM, the following is necessary to make ENERGETICS_SFC_PBL work + !call allocate_forcing_type(G, fluxes, ustar=.true., shelf=.true., & + ! press=.true., water=CS%isthermo, heat=CS%isthermo) call allocate_forcing_type(G, fluxes, ustar=.true., shelf=.true., & - press=.true., water=CS%isthermo, heat=CS%isthermo) + press=.true., water=.true., heat=.true.) ! This is to make it work with layer mode without BULKMIXEDLAYER - !allocate( fluxes%buoy(isd:ied,jsd:jed) ) - !fluxes%buoy(:,:) = 0.0 + allocate( fluxes%buoy(isd:ied,jsd:jed) ) + fluxes%buoy(:,:) = 0.0 else if (is_root_pe()) print *,"allocating fluxes in solo mode" call allocate_forcing_type(G, fluxes, ustar=.true., shelf=.true., press=.true.) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index d2ce432ec5..64907a28d2 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -20,6 +20,7 @@ module ISOMIP_initialization !*********************************************************************** use MOM_ALE_sponge, only : ALE_sponge_CS, set_up_ALE_sponge_field, initialize_ALE_sponge +use MOM_sponge, only : sponge_CS, set_up_sponge_field, initialize_sponge use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type From 80f8b62c6b27e52839256f4ae688f9cebfbfaa7b Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 1 Jul 2016 17:11:32 -0400 Subject: [PATCH 18/46] Flipped geoLat and geoLon to make the domain consistent with the ISOMIP specification. --- src/user/ISOMIP_initialization.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 64907a28d2..f91295d381 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -108,7 +108,7 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) if (is_2D) then do j=js,je ; do i=is,ie ! 2D setup - xtil = G%geoLonT(i,j)*1.0e3/xbar + xtil = G%geoLatT(i,j)*1.0e3/xbar !xtil = 450*1.0e3/xbar bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 !by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & @@ -127,18 +127,18 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) do j=js,je ; do i=is,ie ! 3D setup ! #### TEST ####### - !if (G%geoLonT(i,j)<500.) then + !if (G%geoLatT(i,j)<500.) then ! xtil = 500.*1.0e3/xbar !else - ! xtil = G%geoLonT(i,j)*1.0e3/xbar + ! xtil = G%geoLatT(i,j)*1.0e3/xbar !endif ! ##### TEST ##### - xtil = G%geoLonT(i,j)*1.0e3/xbar + xtil = G%geoLatT(i,j)*1.0e3/xbar bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 - by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & - (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) + by = (dc/(1.+exp(-2.*(G%geoLonT(i,j)*1.0e3- ly/2. - wc)/fc))) + & + (dc/(1.+exp(2.*(G%geoLonT(i,j)*1.0e3- ly/2. + wc)/fc))) D(i,j) = -max(bx+by,-bmax) if (D(i,j) > max_depth) D(i,j) = max_depth @@ -476,10 +476,10 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) ! and mask2dT is 1. do i=is,ie; do j=js,je - if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then + if (G%geoLatT(i,j) >= 790.0 .AND. G%geoLatT(i,j) <= 800.0) then ! 1 / day - dummy1=(G%geoLonT(i,j)-790.0)/(800.0-790.0) + dummy1=(G%geoLatT(i,j)-790.0)/(800.0-790.0) damp = 1.0/TNUDG * max(0.0,dummy1) else ; damp=0.0 From 84c37aa5b336a9673486bd2d8611f59fbe701c51 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 2 Aug 2016 13:20:53 -0400 Subject: [PATCH 19/46] Changed ISOMIP initialization The initial "fit" of temperature or salinity is now done in the ISOMIP initialization routine. This is done in a slighly different way compared to the "fit" option in the MOM initialization file. The relevant parameters are: FIT_SALINITY, DRHO_DS, DRHO_DT, T_REF and S_REF. --- src/initialization/MOM_state_initialization.F90 | 2 +- src/user/ISOMIP_initialization.F90 | 13 ++++++++++--- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 904e2a9af9..c450adb2a5 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1264,7 +1264,7 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, param_file, eqn_of_state, P_Ref units="degC", fail_if_missing=.true.) call get_param(param_file, mod, "S_REF", S_Ref, & "A reference salinity used in initialization.", units="PSU", & - default=35.0) + default=33.8) call get_param(param_file, mod, "FIT_SALINITY", fit_salin, & "If true, accept the prescribed temperature and fit the \n"//& "salinity; otherwise take salinity and fit temperature.", & diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index f41039d6d7..d3dbcfff9f 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -273,7 +273,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, real :: drho_dS(SZK_(G)) ! Derivative of density with salinity in kg m-3 PSU-1. ! real :: rho_guess(SZK_(G)) ! Potential density at T0 & S0 in kg m-3. real :: pres(SZK_(G)) ! An array of the reference pressure in Pa. (zero here) - real :: drho_dT1, drho_dS1 + real :: drho_dT1, drho_dS1, T_Ref, S_Ref is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke pres(:) = 0.0 @@ -319,7 +319,13 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, call get_param(param_file, mod, "DRHO_DT", drho_dT1, & "Partial derivative of density with temperature.", & units="kg m-3 K-1", fail_if_missing=.true.) - + call get_param(param_file, mod, "T_REF", T_Ref, & + "A reference temperature used in initialization.", & + units="degC", fail_if_missing=.true.) + call get_param(param_file, mod, "S_REF", S_Ref, & + "A reference salinity used in initialization.", units="PSU", & + default=35.0) + !write(*,*)'read drho_dS, drho_dT', drho_dS1, drho_dT1 S_range = s_bot - s_sur @@ -330,6 +336,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, i=G%iec; j=G%jec; xi0 = 0.0; do k = 1,nz + !T0(k) = T_Ref; S0(k) = S_Ref xi1 = xi0 + 0.5 * h(i,j,k); S0(k) = S_sur + S_range * xi1; T0(k) = T_sur + T_range * xi1; @@ -370,7 +377,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, enddo endif - do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + do k=1,nz ; do j=js,je ; do i=is,ie T(i,j,k) = T0(k) ; S(i,j,k) = S0(k) enddo ; enddo ; enddo From 826b01431427088506087603d2ce665f93116928 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 2 Aug 2016 16:23:34 -0400 Subject: [PATCH 20/46] Improve remapping in ALE sponge Read parameters (remapping scheme and bnd_extrapolation) from MOM_input so that initialize_remapping is called correctly (i.e., consitent with the way tracers and thickness are remapped in the interior). --- src/parameterizations/vertical/MOM_ALE_sponge.F90 | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index ae2b54f41d..2e074e2ab2 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -84,6 +84,7 @@ subroutine initialize_ALE_sponge(Iresttime, data_h, nz_data, G, param_file, CS) #include "version_variable.h" character(len=40) :: mod = "MOM_sponge" ! This module's name. logical :: use_sponge + logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries integer :: i, j, k, col, total_sponge_cols character(len=10) :: remapScheme if (associated(CS)) then @@ -108,6 +109,13 @@ subroutine initialize_ALE_sponge(Iresttime, data_h, nz_data, G, param_file, CS) " for vertical remapping for all variables.", & default="PLM", do_not_log=.true.) + call get_param(param_file, mod, "BOUNDARY_EXTRAPOLATION", bndExtrapolation, & + "When defined, a proper high-order reconstruction \n"//& + "scheme is used within boundary cells rather \n"// & + "than PCM. E.g., if PPM is used for remapping, a \n" //& + "PPM reconstruction will also be used within boundary cells.", & + default=.false., do_not_log=.true.) + CS%nz = G%ke CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed @@ -149,7 +157,7 @@ subroutine initialize_ALE_sponge(Iresttime, data_h, nz_data, G, param_file, CS) call sum_across_PEs(total_sponge_cols) ! Call the constructor for remapping control structure - call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=.false.) + call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=bndExtrapolation) call log_param(param_file, mod, "!Total sponge columns", total_sponge_cols, & "The total number of columns where sponges are applied.") From 46d8b5d7ac704044f0b068e9f92e9ffb248eed3f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 5 Aug 2016 15:14:49 -0400 Subject: [PATCH 21/46] Added some print statments for debugging purposes --- src/parameterizations/vertical/MOM_ALE_sponge.F90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 2e074e2ab2..1ca67b0018 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -249,9 +249,13 @@ subroutine apply_ALE_sponge(h, dt, G, CS) ! CS%var(m)%p(i,j,k) = I1pdamp * & ! (CS%var(m)%p(i,j,k) + tmp_val1 * damp) enddo - + enddo ! end of c loop - + ! for debugging + !c=CS%num_col + !do m=1,CS%fldno + ! write(*,*)'APPLY SPONGE,m,CS%Ref_h(:,c),h(i,j,:),tmp_val2,tmp_val1',m,CS%Ref_h(:,c),h(i,j,:),tmp_val2,tmp_val1 + !enddo end subroutine apply_ALE_sponge !> GMM: I could not find where sponge_end is being called, but I am keeping From f22c87a66bb8edfc6071e537a7850569ebc6fad4 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 8 Aug 2016 15:33:57 -0400 Subject: [PATCH 22/46] Modified the ice shelf code following the ISOMIP+ setup I've added some of the diagnostics required for the ISOMIP project. I also dded new parameters that can modify the original code: ** the ice shelf can be a perfect insulator (SHELF_INSULATOR = True); ** heat-transfer coeff. can be specified by the user (SHELF_3EQ_GAMMA = True and SHELF_3EQ_GAMMA_T = some value). ** under the ice shelf, ustar is now computed following the ISOMIP setup ustar^2 = Cd*(u^2 + u_tidal^2) --- src/ice_shelf/MOM_ice_shelf.F90 | 102 ++++++++++++++++++++++++-------- 1 file changed, 77 insertions(+), 25 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 1136b32941..228cf75fc0 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -255,7 +255,8 @@ module MOM_ice_shelf real :: kd_molec_salt ! The molecular diffusivity of salt, in m2 s-1. real :: kd_molec_temp ! The molecular diffusivity of heat, in m2 s-1. real :: Lat_fusion ! The latent heat of fusion, in J kg-1. - + real :: Gamma_T_3EQ ! Nondimensional heat-transfer coefficient, used in the 3Eq. formulation + ! This number should be specified by the user. real :: col_thick_melt_threshold ! if the mixed layer is below this threshold, melt rate ! is not calculated for the cell @@ -334,13 +335,18 @@ module MOM_ice_shelf ! the underlying ocean. logical :: threeeq ! If true, the 3 equation consistency equations are ! used to calculate the flux at the ocean-ice interface. + logical :: insulator ! If true, ice shelf is a perfect insulator + logical :: const_gamma ! If true, gamma_T is specified by the user. + integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, & id_tfreeze = -1, id_tfl_shelf = -1, & - id_u_shelf = -1, id_v_shelf = -1, id_h_shelf = -1, id_h_mask = -1, & - id_u_mask = -1, id_v_mask = -1, id_t_shelf = -1, id_t_mask = -1, & - id_surf_elev = -1, id_bathym = -1, id_float_frac = -1, id_col_thick = -1, & - id_area_shelf_h = -1, id_OD_av = -1, id_float_frac_rt = -1, id_ustar_shelf = -1, & - id_shelf_mass = -1 + id_thermal_driving = -1, id_haline_driving = -1, & + id_u_ml = -1, id_v_ml = -1, & + id_u_shelf = -1, id_v_shelf = -1, id_h_shelf = -1, id_h_mask = -1, & + id_u_mask = -1, id_v_mask = -1, id_t_shelf = -1, id_t_mask = -1, & + id_surf_elev = -1, id_bathym = -1, id_float_frac = -1, id_col_thick = -1, & + id_area_shelf_h = -1, id_OD_av = -1, id_float_frac_rt = -1,& + id_ustar_shelf = -1, id_shelf_mass = -1 ! ids for outputting intermediate thickness in advection subroutine (debugging) !integer :: id_h_after_uflux = -1, id_h_after_vflux = -1, id_h_after_adv = -1 @@ -495,8 +501,9 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) iDens = 1.0/CS%density_ocean_avg - if (CS%shelf_mass_is_dynamic .and. CS%override_shelf_movement) call update_shelf_mass(CS, Time) + if (CS%shelf_mass_is_dynamic .and. CS%override_shelf_movement) call update_shelf_mass(CS, Time) + do j=js,je ! Find the pressure at the ice-ocean interface, averaged only over the ! part of the cell covered by ice shelf. @@ -525,8 +532,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) u_at_h = state%u(i,j) v_at_h = state%v(i,j) - fluxes%ustar_shelf(i,j)= sqrt(CS%cdrag)*((u_at_h**2.0 + v_at_h**2.0)**0.5 +& - CS%utide(i,j)) + fluxes%ustar_shelf(i,j)= sqrt(CS%cdrag)*sqrt((u_at_h**2.0 + v_at_h**2.0)**2 +& + CS%utide(i,j)**2) ustar_h = MAX(CS%ustar_bg, fluxes%ustar_shelf(i,j)) fluxes%ustar_shelf(i,j) = ustar_h @@ -563,10 +570,18 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) ! First, determine the buoyancy flux assuming no effects of stability ! on the turbulence. Following H & J '99, this limit also applies ! when the buoyancy flux is destabilizing. + + if (CS%const_gamma) then ! if using a constant gamma_T + !I_Gam_T = 1.0 / CS%Gamma_T_3EQ + !I_Gam_S = 1.0 / (CS%Gamma_T_3EQ/35.) + I_Gam_T = CS%Gamma_T_3EQ + I_Gam_S = CS%Gamma_T_3EQ/35. + else + Gam_turb = I_VK * (ln_neut + (0.5 * I_ZETA_N - 1.0)) + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + endif - Gam_turb = I_VK * (ln_neut + (0.5 * I_ZETA_N - 1.0)) - I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) - I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) wT_flux = dT_ustar * I_Gam_T wB_flux = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * wT_flux @@ -592,8 +607,16 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) dG_dwB = I_VK * (0.5 * I_ZETA_N) * dIns_dwB endif - I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) - I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + if (CS%const_gamma) then ! if using a constant gamma_T + !I_Gam_T = 1.0 / CS%Gamma_T_3EQ + I_Gam_T = CS%Gamma_T_3EQ + !I_Gam_S = 1.0 / (CS%Gamma_T_3EQ/35.) + I_Gam_S = CS%Gamma_T_3EQ/35. + else + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + endif + wT_flux = dT_ustar * I_Gam_T wB_flux_new = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * wT_flux @@ -626,21 +649,26 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) CS%lprec(i,j) = I_LF * CS%t_flux(i,j) CS%tflux_shelf(i,j) = 0.0 else - ! With melting, from H&J 1999, eqs (31) & (26)... - ! Q_ice ~= cp_ice * (CS%Temp_Ice-T_freeze) * lprec - ! RhoLF*lprec = Q_ice + CS%t_flux(i,j) - ! lprec = (CS%t_flux(i,j)) / (LF + cp_ice * (T_freeze-CS%Temp_Ice)) - CS%lprec(i,j) = CS%t_flux(i,j) / & - (LF + CS%CP_Ice * (CS%Tfreeze(i,j) - CS%Temp_Ice)) - - CS%tflux_shelf(i,j) = CS%t_flux(i,j) - LF*CS%lprec(i,j) + if (CS%insulator) then + !no conduction/perfect insulator + CS%tflux_shelf(i,j) = 0.0 + CS%lprec(i,j) = I_LF * (- CS%tflux_shelf(i,j) + CS%t_flux(i,j)) + + else + ! With melting, from H&J 1999, eqs (31) & (26)... + ! Q_ice ~= cp_ice * (CS%Temp_Ice-T_freeze) * lprec + ! RhoLF*lprec = Q_ice + CS%t_flux(i,j) + ! lprec = (CS%t_flux(i,j)) / (LF + cp_ice * (T_freeze-CS%Temp_Ice)) + CS%lprec(i,j) = CS%t_flux(i,j) / & + (LF + CS%CP_Ice * (CS%Tfreeze(i,j) - CS%Temp_Ice)) + + CS%tflux_shelf(i,j) = CS%t_flux(i,j) - LF*CS%lprec(i,j) + endif + endif !other options: dTi/dz linear through shelf ! dTi_dz = (CS%Temp_Ice - CS%tfreeze(i,j))/G%draft(i,j) ! CS%tflux_shelf(i,j) = - Rho_Ice * CS%CP_Ice * KTI * dTi_dz - !no conduction/perfect insulator - ! CS%tflux_shelf(i,j) = 0.0 - ! CS%lprec(i,j) = I_LF * (- CS%tflux_shelf(i,j) + CS%t_flux(i,j)) mass_exch = CS%exch_vel_s(i,j) * CS%Rho0 Sbdry_it = (state%sss(i,j) * mass_exch + CS%Salin_ice * CS%lprec(i,j)) / & @@ -688,6 +716,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) ! melt in m/year fluxes%iceshelf_melt = CS%lprec * (86400.0*365.0/CS%density_ice) + !CS%lprec(:,:) = 5.822E-6 ! testing if (CS%DEBUG) then call hchksum (CS%h_shelf, "melt rate", G%HI, haloshift=0) @@ -750,6 +779,10 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) if (CS%id_area_shelf_h > 0) call post_data(CS%id_area_shelf_h, CS%area_shelf_h, CS%diag) if (CS%id_ustar_shelf > 0) call post_data(CS%id_ustar_shelf, fluxes%ustar_shelf, CS%diag) if (CS%id_melt > 0) call post_data(CS%id_melt, (CS%lprec) * (86400.0*365.0/CS%density_ice), CS%diag) + if (CS%id_thermal_driving > 0) call post_data(CS%id_thermal_driving, (state%sst-CS%tfreeze), CS%diag) + if (CS%id_haline_driving > 0) call post_data(CS%id_haline_driving, (state%sss-Sbdry), CS%diag) + if (CS%id_u_ml > 0) call post_data(CS%id_u_ml,state%u,CS%diag) + if (CS%id_v_ml > 0) call post_data(CS%id_v_ml,state%v,CS%diag) if (CS%id_tfreeze > 0) call post_data(CS%id_tfreeze, CS%tfreeze, CS%diag) if (CS%id_tfl_shelf > 0) call post_data(CS%id_tfl_shelf, CS%tflux_shelf, CS%diag) if (CS%id_exch_vel_t > 0) call post_data(CS%id_exch_vel_t, CS%exch_vel_t, CS%diag) @@ -1160,6 +1193,17 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti "If true, use the three equation expression of \n"//& "consistency to calculate the fluxes at the ice-ocean \n"//& "interface.", default=.true.) + call get_param(param_file, mod, "SHELF_INSULATOR", CS%insulator, & + "If true, the ice shelf is a perfect insulatior \n"//& + "(no conduction).", default=.false.) + call get_param(param_file, mod, "SHELF_3EQ_GAMMA", CS%const_gamma, & + "If true, user specifies a constant nondimensional heat-transfer coefficient \n"//& + "(GAMMA_T_3EQ), from which the salt-transfer coefficient is then computed \n"//& + " as GAMMA_T_3EQ/35. This is used with SHELF_THREE_EQN.", default=.false.) + if (CS%const_gamma) call get_param(param_file, mod, "SHELF_3EQ_GAMMA_T", CS%Gamma_T_3EQ, & + "Nondimensional heat-transfer coefficient.",default=2.2E-2, & + units="nondim.", fail_if_missing=.true.) + if (.not.CS%threeeq) & call get_param(param_file, mod, "SHELF_2EQ_GAMMA_T", CS%gamma_t, & "If SHELF_THREE_EQN is false, this the fixed turbulent \n"//& @@ -1688,6 +1732,14 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti 'mass of shelf', 'kg/m^2') CS%id_melt = register_diag_field('ocean_model', 'melt', CS%diag%axesT1, CS%Time, & 'Ice Shelf Melt Rate', 'meter year-1') + CS%id_thermal_driving = register_diag_field('ocean_model', 'thermal_driving', CS%diag%axesT1, CS%Time, & + 'pot. temp. in the boundary layer minus freezing pot. temp. at the ice-ocean interface.', 'Celsius') + CS%id_haline_driving = register_diag_field('ocean_model', 'haline_driving', CS%diag%axesT1, CS%Time, & + 'salinity in the boundary layer minus salinity at the ice-ocean interface.', 'PPT') + CS%id_u_ml = register_diag_field('ocean_model', 'u_ml', CS%diag%axesT1, CS%Time, & + 'Eastward vel. in the boundary layer (used to compute ustar)', 'meter second-1') + CS%id_v_ml = register_diag_field('ocean_model', 'v_ml', CS%diag%axesT1, CS%Time, & + 'Northward vel. in the boundary layer (used to compute ustar)', 'meter second-1') CS%id_exch_vel_s = register_diag_field('ocean_model', 'exch_vel_s', CS%diag%axesT1, CS%Time, & 'Sub-shelf salinity exchange velocity', 'meter second-1') CS%id_exch_vel_t = register_diag_field('ocean_model', 'exch_vel_t', CS%diag%axesT1, CS%Time, & From 98824eea5cd54efe1494180923df3fd3404cf60c Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 8 Aug 2016 15:44:14 -0400 Subject: [PATCH 23/46] Added some write statements in ISOMIP_initialization for debugging --- src/user/ISOMIP_initialization.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index d3dbcfff9f..421596c95c 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -582,7 +582,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) !i=G%iec; j=G%jec !do k = 1,nz ! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state) - !write(*,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) + ! write(*,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) !enddo ! Now register all of the fields which are damped in the sponge. ! From e7894ad29805c0dbfecfab8670bcf87c410a58d3 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 15 Aug 2016 13:51:18 -0400 Subject: [PATCH 24/46] Added new diagnostics for ISOMIP --- src/ice_shelf/MOM_ice_shelf.F90 | 33 ++++++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 9 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 228cf75fc0..2a0d79de3d 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -346,7 +346,7 @@ module MOM_ice_shelf id_u_mask = -1, id_v_mask = -1, id_t_shelf = -1, id_t_mask = -1, & id_surf_elev = -1, id_bathym = -1, id_float_frac = -1, id_col_thick = -1, & id_area_shelf_h = -1, id_OD_av = -1, id_float_frac_rt = -1,& - id_ustar_shelf = -1, id_shelf_mass = -1 + id_ustar_shelf = -1, id_shelf_mass = -1, id_mass_flux = -1 ! ids for outputting intermediate thickness in advection subroutine (debugging) !integer :: id_h_after_uflux = -1, id_h_after_vflux = -1, id_h_after_adv = -1 @@ -433,7 +433,12 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) ! with temperature, in units of kg m-3 K-1. dR0_dS, & ! Partial derivative of the mixed layer density ! with salinity, in units of kg m-3 psu-1. - p_int ! The pressure at the ice-ocean interface, in Pa. + p_int ! The pressure at the ice-ocean interface, in Pa. + + real, dimension(:,:), allocatable :: mass_flux ! total mass flux of freshwater across + ! ice-ocean interface, positive for melting + ! and negative for freezing. This is computed + ! as part of the ISOMIP diagnostics. real, parameter :: VK = 0.40 ! Von Karman's constant - dimensionless real :: ZETA_N = 0.052 ! The fraction of the boundary layer over which the ! viscosity is linearly increasing. (Was 1/8. Why?) @@ -477,7 +482,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) type(ocean_grid_type), pointer :: G real, parameter :: c2_3 = 2.0/3.0 integer :: i, j, is, ie, js, je, ied, jed, it1, it3, iters_vel_solve - + real, parameter :: rho_fw = 1000.0 ! fresh water density if (.not. associated(CS)) call MOM_error(FATAL, "shelf_calc_flux: "// & "initialize_ice_shelf must be called before shelf_calc_flux.") call cpu_clock_begin(id_clock_shelf) @@ -714,10 +719,16 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) enddo ! i-loop enddo ! j-loop - ! melt in m/year - fluxes%iceshelf_melt = CS%lprec * (86400.0*365.0/CS%density_ice) - !CS%lprec(:,:) = 5.822E-6 ! testing - + ! CS%lprec = precipitating liquid water into the ocean ( kg/(m^2 s) ) + ! melt in m/year + !fluxes%iceshelf_melt = CS%lprec * (86400.0*365.0/CS%density_ice) + ! Use rho_fw for ISOMIP experiments + fluxes%iceshelf_melt = CS%lprec * (86400.0*365.0/rho_fw) + ! mass flux (kg/s), part of ISOMIP diags. + ! first, allocate mass_flux. This is probably not the best way to do this. + !im,jm=SHAPE(CS%lprec) + ALLOCATE ( mass_flux(G%ied,G%jed) ) + mass_flux = (CS%lprec) * CS%area_shelf_h if (CS%DEBUG) then call hchksum (CS%h_shelf, "melt rate", G%HI, haloshift=0) endif @@ -778,9 +789,12 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) if (CS%id_shelf_mass > 0) call post_data(CS%id_shelf_mass, CS%mass_shelf, CS%diag) if (CS%id_area_shelf_h > 0) call post_data(CS%id_area_shelf_h, CS%area_shelf_h, CS%diag) if (CS%id_ustar_shelf > 0) call post_data(CS%id_ustar_shelf, fluxes%ustar_shelf, CS%diag) - if (CS%id_melt > 0) call post_data(CS%id_melt, (CS%lprec) * (86400.0*365.0/CS%density_ice), CS%diag) + ! Use freshwater density fotr ISOMIP exps. + if (CS%id_melt > 0) call post_data(CS%id_melt, (CS%lprec) * (86400.0*365.0/rho_fw), CS%diag) + !if (CS%id_melt > 0) call post_data(CS%id_melt, (CS%lprec) * (86400.0*365.0/CS%density_ice), CS%diag) if (CS%id_thermal_driving > 0) call post_data(CS%id_thermal_driving, (state%sst-CS%tfreeze), CS%diag) if (CS%id_haline_driving > 0) call post_data(CS%id_haline_driving, (state%sss-Sbdry), CS%diag) + if (CS%id_mass_flux > 0) call post_data(CS%id_mass_flux, mass_flux, CS%diag) if (CS%id_u_ml > 0) call post_data(CS%id_u_ml,state%u,CS%diag) if (CS%id_v_ml > 0) call post_data(CS%id_v_ml,state%v,CS%diag) if (CS%id_tfreeze > 0) call post_data(CS%id_tfreeze, CS%tfreeze, CS%diag) @@ -1109,7 +1123,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti logical :: solo_ice_sheet, read_TideAmp character(len=240) :: Tideamp_file real :: utide - if (associated(CS)) then call MOM_error(FATAL, "MOM_ice_shelf.F90, initialize_ice_shelf: "// & "called with an associated control structure.") @@ -1730,6 +1743,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti 'Ice Shelf Area in cell', 'meter-2') CS%id_shelf_mass = register_diag_field('ocean_model', 'shelf_mass', CS%diag%axesT1, CS%Time, & 'mass of shelf', 'kg/m^2') + CS%id_mass_flux = register_diag_field('ocean_model', 'mass_flux', CS%diag%axesT1,& + CS%Time,'Total mass flux of freshwater across the ice-ocean interface.', 'kg/s') CS%id_melt = register_diag_field('ocean_model', 'melt', CS%diag%axesT1, CS%Time, & 'Ice Shelf Melt Rate', 'meter year-1') CS%id_thermal_driving = register_diag_field('ocean_model', 'thermal_driving', CS%diag%axesT1, CS%Time, & From d68270da8d6e15843cbb632d66ea27c72660f4c4 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 16 Aug 2016 20:55:54 -0400 Subject: [PATCH 25/46] Changed ISOMIP initialization for layer mode Fixed a bug in the initial layer "fit" of temperature or salinity done in the ISOMIP initialization routine. --- src/user/ISOMIP_initialization.F90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 421596c95c..1bb1fa4b85 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -334,7 +334,8 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, S_range = S_range / G%max_depth ! Convert S_range into dS/dz T_range = T_range / G%max_depth ! Convert T_range into dT/dz - i=G%iec; j=G%jec; xi0 = 0.0; + do j=js,je ; do i=is,ie + xi0 = 0.0; do k = 1,nz !T0(k) = T_Ref; S0(k) = S_Ref xi1 = xi0 + 0.5 * h(i,j,k); @@ -377,9 +378,11 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, enddo endif - do k=1,nz ; do j=js,je ; do i=is,ie + do k=1,nz T(i,j,k) = T0(k) ; S(i,j,k) = S0(k) - enddo ; enddo ; enddo + enddo + + enddo ; enddo case default call MOM_error(FATAL,"isomip_initialize: "// & From e2eb3b3adc0fe6031cf5badd26d391f04a576322 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 18 Aug 2016 12:49:40 -0400 Subject: [PATCH 26/46] Added parameter MELTING_CUTOFF_DEPTH This sets the depth above which the melt is set to zero (required for the ISOMIP test case). I've used pressure (mass of ice shelf * RHO_0 * g) to find out whether lprec needs to be set to zero. Since RHO_0 is used, this is an approximation which works fine in this case. --- src/ice_shelf/MOM_ice_shelf.F90 | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 2a0d79de3d..8d3ed9f109 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -337,7 +337,7 @@ module MOM_ice_shelf ! used to calculate the flux at the ocean-ice interface. logical :: insulator ! If true, ice shelf is a perfect insulator logical :: const_gamma ! If true, gamma_T is specified by the user. - + real :: cutoff_depth ! depth above which melt is set to zero (>= 0). integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, & id_tfreeze = -1, id_tfl_shelf = -1, & id_thermal_driving = -1, id_haline_driving = -1, & @@ -507,7 +507,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) iDens = 1.0/CS%density_ocean_avg - if (CS%shelf_mass_is_dynamic .and. CS%override_shelf_movement) call update_shelf_mass(CS, Time) + if (CS%shelf_mass_is_dynamic .and. CS%override_shelf_movement) & + call update_shelf_mass(CS, Time) do j=js,je ! Find the pressure at the ice-ocean interface, averaged only over the @@ -525,7 +526,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) ! DNG - to allow this everywhere Hml>0.0 allows for melting under grounded cells ! propose instead to allow where Hml > [some threshold] - if ((iDens*state%ocean_mass(i,j) > CS%col_thick_melt_threshold) .and. (CS%area_shelf_h(i,j) > 0.0) .and. & + if ((iDens*state%ocean_mass(i,j) > CS%col_thick_melt_threshold) .and. & + (CS%area_shelf_h(i,j) > 0.0) .and. & (CS%isthermo) .and. (state%Hml(i,j) > 0.0) ) then if (CS%threeeq) then @@ -716,6 +718,11 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) else !not shelf CS%t_flux(i,j) = 0.0 endif + + ! GM: set lprec (therefore melting) to zero above a cutoff pressure + ! (CS%Rho0*CS%cutoff_depth*CS%g_Earth) this is needed for the isomip test case + if (p_int(i) < CS%Rho0*CS%cutoff_depth*CS%g_Earth) CS%lprec(i,j) = 0.0 + enddo ! i-loop enddo ! j-loop @@ -1209,6 +1216,12 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti call get_param(param_file, mod, "SHELF_INSULATOR", CS%insulator, & "If true, the ice shelf is a perfect insulatior \n"//& "(no conduction).", default=.false.) + call get_param(param_file, mod, "MELTING_CUTOFF_DEPTH", CS%cutoff_depth, & + "Depth above which the melt is set to zero (it must be >= 0) \n"//& + "Default value won't affect the solution.", default=0.0) + if (CS%cutoff_depth < 0.) & + call MOM_error(WARNING,"Initialize_ice_shelf: MELTING_CUTOFF_DEPTH must be >= 0.") + call get_param(param_file, mod, "SHELF_3EQ_GAMMA", CS%const_gamma, & "If true, user specifies a constant nondimensional heat-transfer coefficient \n"//& "(GAMMA_T_3EQ), from which the salt-transfer coefficient is then computed \n"//& @@ -1473,6 +1486,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti call register_restart_field(CS%mass_shelf, vd, .true., CS%restart_CSp) vd = var_desc("shelf_area","m2","Ice shelf area in cell",z_grid='1') call register_restart_field(CS%area_shelf_h, vd, .true., CS%restart_CSp) + vd = var_desc("h_shelf","m","ice sheet/shelf thickness",z_grid='1') + call register_restart_field(CS%h_shelf, vd, .true., CS%restart_CSp) if (CS%shelf_mass_is_dynamic .and. .not.CS%override_shelf_movement) then ! additional restarts for ice shelf state @@ -1480,8 +1495,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti call register_restart_field(CS%u_shelf, vd, .true., CS%restart_CSp) vd = var_desc("v_shelf","m s-1","ice sheet/shelf velocity",'q',z_grid='1') call register_restart_field(CS%v_shelf, vd, .true., CS%restart_CSp) - vd = var_desc("h_shelf","m","ice sheet/shelf thickness",z_grid='1') - call register_restart_field(CS%h_shelf, vd, .true., CS%restart_CSp) + !vd = var_desc("h_shelf","m","ice sheet/shelf thickness",z_grid='1') + !call register_restart_field(CS%h_shelf, vd, .true., CS%restart_CSp) vd = var_desc("h_mask","none","ice sheet/shelf thickness mask",z_grid='1') call register_restart_field(CS%hmask, vd, .true., CS%restart_CSp) @@ -1764,7 +1779,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti CS%id_tfl_shelf = register_diag_field('ocean_model', 'tflux_shelf', CS%diag%axesT1, CS%Time, & 'Heat conduction into ice shelf', 'Watts meter-2') CS%id_ustar_shelf = register_diag_field('ocean_model', 'ustar_shelf', CS%diag%axesT1, CS%Time, & - 'Fric vel under shelf', 'm/s?') + 'Fric vel under shelf', 'm/s') + if (CS%shelf_mass_is_dynamic .and. .not.CS%override_shelf_movement) then CS%id_u_shelf = register_diag_field('ocean_model','u_shelf',CS%diag%axesB1,CS%Time, & 'x-velocity of ice', 'm year') @@ -1774,8 +1790,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti 'mask for u-nodes', 'none') CS%id_v_mask = register_diag_field('ocean_model','v_mask',CS%diag%axesB1,CS%Time, & 'mask for v-nodes', 'none') - CS%id_h_shelf = register_diag_field('ocean_model','h_shelf',CS%diag%axesT1,CS%Time, & - 'land ice thickness', 'm') CS%id_h_mask = register_diag_field('ocean_model','h_mask',CS%diag%axesT1,CS%Time, & 'ice shelf thickness', 'none') CS%id_surf_elev = register_diag_field('ocean_model','ice_surf',CS%diag%axesT1,CS%Time, & From 081b19d96ef6a764347a14a62f9fdc6b6b108044 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 25 Aug 2016 13:50:01 -0400 Subject: [PATCH 27/46] Added a new paramter in the ice shelf code and modified isomip_tracer Added a new parameter (MELTING_CUTOFF_DEPTH) in MOM_ice_shelf.F90 that controls the above which the melt is set to zero. Modified ISOMIP_tracer.F90 so that only one tracer is now injected under the ice shelf. There was another tracer being injected in the sponge layer (code for that is not commented), but that did not wrok in layer mode (it needs to be modified to allow sponge in addition to sponge_ale). --- src/ice_shelf/MOM_ice_shelf.F90 | 26 +++++++++---- src/tracer/ISOMIP_tracer.F90 | 67 +++++++++++++++++---------------- 2 files changed, 53 insertions(+), 40 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 8d3ed9f109..9376cf8e60 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -489,7 +489,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) G => CS%grid is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; ied = G%ied ; jed = G%jed - + Sbdry = 0.0 ! define Sbdry to avoid Run-Time Check Failure, when melt is not computed. I_ZETA_N = 1.0 / ZETA_N LF = CS%Lat_fusion I_RhoLF = 1.0/(CS%Rho0*LF) @@ -699,7 +699,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) Sbdry = Sbdry_it endif - Sbdry = Sbdry_it + Sbdry = Sbdry_it enddo !it1 ! Check for non-convergence and/or non-boundedness? @@ -714,15 +714,26 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) CS%t_flux(i,j) = RhoCp * CS%exch_vel_t(i,j) * (state%sst(i,j) - CS%tfreeze(i,j)) CS%tflux_shelf(i,j) = 0.0 CS%lprec(i,j) = I_LF * CS%t_flux(i,j) + Sbdry = 0.0 endif else !not shelf CS%t_flux(i,j) = 0.0 endif - ! GM: set lprec (therefore melting) to zero above a cutoff pressure - ! (CS%Rho0*CS%cutoff_depth*CS%g_Earth) this is needed for the isomip test case - if (p_int(i) < CS%Rho0*CS%cutoff_depth*CS%g_Earth) CS%lprec(i,j) = 0.0 + enddo ! i-loop + enddo ! j-loop + do j=js,je + do i=is,ie + ! GM: set lprec (therefore melting) to zero above a cutoff pressure + ! (CS%Rho0*CS%cutoff_depth*CS%g_Earth) this is needed for the isomip + ! test case. + + if ((CS%g_Earth * CS%mass_shelf(i,j)) < CS%Rho0*CS%cutoff_depth* & + CS%g_Earth) CS%lprec(i,j) = 0.0 + ! avoid negative fluxes + ! CS%lprec(i,j) = MAX(CS%lprec(i,j),0.0) + enddo ! i-loop enddo ! j-loop @@ -800,7 +811,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) if (CS%id_melt > 0) call post_data(CS%id_melt, (CS%lprec) * (86400.0*365.0/rho_fw), CS%diag) !if (CS%id_melt > 0) call post_data(CS%id_melt, (CS%lprec) * (86400.0*365.0/CS%density_ice), CS%diag) if (CS%id_thermal_driving > 0) call post_data(CS%id_thermal_driving, (state%sst-CS%tfreeze), CS%diag) - if (CS%id_haline_driving > 0) call post_data(CS%id_haline_driving, (state%sss-Sbdry), CS%diag) + if (CS%id_haline_driving > 0) call post_data(CS%id_haline_driving, (state%sss - Sbdry), CS%diag) if (CS%id_mass_flux > 0) call post_data(CS%id_mass_flux, mass_flux, CS%diag) if (CS%id_u_ml > 0) call post_data(CS%id_u_ml,state%u,CS%diag) if (CS%id_v_ml > 0) call post_data(CS%id_v_ml,state%v,CS%diag) @@ -1413,6 +1424,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti allocate( CS%tflux_shelf(isd:ied,jsd:jed) ) ; CS%tflux_shelf(:,:) = 0.0 allocate( CS%tfreeze(isd:ied,jsd:jed) ) ; CS%tfreeze(:,:) = 0.0 +! allocate( CS%Sbdry(isd:ied,jsd:jed) ) ; CS%Sbdry(:,:) = 0.0 allocate( CS%exch_vel_s(isd:ied,jsd:jed) ) ; CS%exch_vel_s(:,:) = 0.0 allocate( CS%exch_vel_t(isd:ied,jsd:jed) ) ; CS%exch_vel_t(:,:) = 0.0 @@ -5614,7 +5626,7 @@ subroutine ice_shelf_end(CS) deallocate(CS%t_flux) ; deallocate(CS%lprec) deallocate(CS%salt_flux) - deallocate(CS%tflux_shelf) ; deallocate(CS%tfreeze) + deallocate(CS%tflux_shelf) ; deallocate(CS%tfreeze);! deallocate(CS%Sbdry); deallocate(CS%exch_vel_t) ; deallocate(CS%exch_vel_s) deallocate(CS%h_shelf) ; deallocate(CS%hmask) diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index 6b826f40d5..027e06389d 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -46,7 +46,7 @@ module ISOMIP_tracer public ISOMIP_tracer_column_physics, ISOMIP_tracer_surface_state, ISOMIP_tracer_end !< ntr is the number of tracers in this module. -integer, parameter :: ntr = 2 +integer, parameter :: ntr = 1 type p3d real, dimension(:,:,:), pointer :: p => NULL() @@ -225,40 +225,41 @@ subroutine initialize_ISOMIP_tracer(restart, day, G, GV, h, diag, OBC, CS, & else do m=1,NTR do k=1,nz ; do j=js,je ; do i=is,ie - CS%tr(i,j,k,m) = 1.0e-20 ! This could just as well be 0. + CS%tr(i,j,k,m) = 0.0 enddo ; enddo ; enddo enddo endif endif ! restart - if ( CS%use_sponge ) then -! If sponges are used, this example damps tracers in sponges in the -! northern half of the domain to 1 and tracers in the southern half -! to 0. For any tracers that are not damped in the sponge, the call -! to set_up_sponge_field can simply be omitted. - if (.not.associated(ALE_sponge_CSp)) & - call MOM_error(FATAL, "ISOMIP_initialize_tracer: "// & - "The pointer to ALEsponge_CSp must be associated if SPONGE is defined.") - - allocate(temp(G%isd:G%ied,G%jsd:G%jed,nz)) +! the following does not work in layer mode yet +!! if ( CS%use_sponge ) then + ! If sponges are used, this example damps tracers in sponges in the + ! northern half of the domain to 1 and tracers in the southern half + ! to 0. For any tracers that are not damped in the sponge, the call + ! to set_up_sponge_field can simply be omitted. +! if (.not.associated(ALE_sponge_CSp)) & +! call MOM_error(FATAL, "ISOMIP_initialize_tracer: "// & +! "The pointer to ALEsponge_CSp must be associated if SPONGE is defined.") + +! allocate(temp(G%isd:G%ied,G%jsd:G%jed,nz)) - do j=js,je ; do i=is,ie - if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then - temp(i,j,:) = 1.0 - else - temp(i,j,:) = 0.0 - endif - enddo ; enddo - -! do m=1,NTR - do m=1,1 +! do j=js,je ; do i=is,ie +! if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then +! temp(i,j,:) = 1.0 +! else +! temp(i,j,:) = 0.0 +! endif +! enddo ; enddo + + ! do m=1,NTR +! do m=1,1 ! This is needed to force the compiler not to do a copy in the sponge ! calls. Curses on the designers and implementers of Fortran90. - tr_ptr => CS%tr(:,:,:,m) - call set_up_ALE_sponge_field(temp, G, tr_ptr, ALE_sponge_CSp) - enddo - deallocate(temp) - endif +! tr_ptr => CS%tr(:,:,:,m) +! call set_up_ALE_sponge_field(temp, G, tr_ptr, ALE_sponge_CSp) +! enddo +! deallocate(temp) +! endif ! This needs to be changed if the units of tracer are changed above. if (GV%Boussinesq) then ; flux_units = "kg kg-1 m3 s-1" @@ -329,6 +330,7 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G ! The arguments to this subroutine are redundant in that ! h_new[k] = h_old[k] + ea[k] - eb[k-1] + eb[k] - ea[k+1] + real :: mmax real :: b1(SZI_(G)) ! b1 and c1 are variables used by the real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver. integer :: i, j, k, is, ie, js, je, nz, m @@ -337,14 +339,13 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G if (.not.associated(CS)) return - ! dye melt water (m=2) - ! converting melt (m/year) to (m/s) - do m=2,NTR + ! max. melt + mmax = MAXVAL(fluxes%iceshelf_melt) + ! dye melt water (m=1), dye = 1 if melt=max(melt) + do m=1,NTR do j=js,je ; do i=is,ie if (fluxes%iceshelf_melt(i,j) > 0.0) then -! CS%tr(i,j,1,m) = (dt*fluxes%iceshelf_melt(i,j)/(365.0 * 86400.0) & -! + h_old(i,j,1)*CS%tr(i,j,1,m))/h_new(i,j,1) - CS%tr(i,j,1:3,m) = 1.0 + CS%tr(i,j,1,m) = fluxes%iceshelf_melt(i,j)/mmax endif enddo ; enddo enddo From b6a06827d9085143e35961c6ffb803174f9d6092 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 25 Aug 2016 17:14:22 -0400 Subject: [PATCH 28/46] Compute velocities (and add to state) in the ML when *not* using BULKMIXLAYER When BULKMIXLAYER is not defined, both u and v are averaged over the uppermost depth_ml fluid. Before, just tracers were being averaged over depth_ml. This change is needed to correctly compute u and v in the ML when layer thicknesses are vanished (i.e., under an ice shelf). --- src/core/MOM.F90 | 95 ++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 84 insertions(+), 11 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c466b8f620..eed2fedd1e 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1533,7 +1533,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call get_param(param_file, "MOM", "HMIX_SFC_PROP", CS%Hmix, & "If BULKMIXEDLAYER is false, HMIX_SFC_PROP is the depth \n"//& "over which to average to find surface properties like \n"//& - "SST and SSS or density (but not surface velocities).", & + "SST, SSS, density are surface velocities.", & units="m", default=1.0) endif call get_param(param_file, "MOM", "MIN_Z_DIAG_INTERVAL", Z_diag_int, & @@ -2855,22 +2855,28 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, GV, CS, p_atm) real, optional, pointer, dimension(:,:) :: p_atm !< atmospheric pressure (Pascal) ! local - real :: depth(SZI_(G)) ! distance from the surface (meter) - real :: depth_ml ! depth over which to average to - ! determine mixed layer properties (meter) - real :: dh ! thickness of a layer within mixed layer (meter) - real :: mass ! mass per unit area of a layer (kg/m2) + real :: hu(SZIB_(G),SZJ_(G),SZK_(G))! layer thickness (m or kg/m2) at u points + real :: hv(SZI_(G),SZJB_(G),SZK_(G))! layer thickness (m or kg/m2) at v points + real :: depth(SZI_(G)) ! distance from the surface (meter) + real :: depth_ml ! depth over which to average to + ! determine mixed layer properties (meter) + real :: dh ! thickness of a layer within mixed layer (meter) + real :: mass ! mass per unit area of a layer (kg/m2) real :: IgR0 integer :: i, j, k, is, ie, js, je, nz, numberOfErrors - integer :: isd, ied, jsd, jed + integer :: isd, ied, jsd, jed + integer :: iscB, iecB, jscB, jecB, isdB, iedB, jsdB, jedB logical :: localError character(240) :: msg call callTree_enter("calculate_surface_state(), MOM.F90") is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + iscB = G%iscB ; iecB = G%iecB; jscB = G%jscB ; jecB = G%jecB + isdB = G%isdB ; iedB = G%iedB; jsdB = G%jsdB ; jedB = G%jedB + !write(*,*)'iscB,iecB,jscB,jecB', iscB,iecB,jscB,jecB state%sea_lev => ssh if (present(p_atm)) then ; if (ASSOCIATED(p_atm)) then @@ -2883,6 +2889,8 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, GV, CS, p_atm) if (CS%bulkmixedlayer) then state%SST => CS%tv%T(:,:,1) state%SSS => CS%tv%S(:,:,1) + state%u => u(:,:,1) + state%v => v(:,:,1) nullify(state%sfc_density) if (associated(CS%tv%Hml)) state%Hml => CS%tv%Hml else @@ -2900,10 +2908,17 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, GV, CS, p_atm) endif nullify(state%SST) ; nullify(state%SSS) endif + if (.not.associated(state%u)) then + allocate(state%u(isdB:iedB,jsd:jed)) ; state%u(:,:) = 0.0 + endif + if (.not.associated(state%v)) then + allocate(state%v(isd:ied,jsdB:jedB)) ; state%v(:,:) = 0.0 + endif + if (.not.associated(state%Hml)) allocate(state%Hml(isd:ied,jsd:jed)) depth_ml = CS%Hmix - ! Determine the mean properties of the uppermost depth_ml fluid. + ! Determine the mean tracer properties of the uppermost depth_ml fluid. !$OMP parallel do default(none) shared(is,ie,js,je,nz,CS,state,h,depth_ml,GV) private(depth,dh) do j=js,je do i=is,ie @@ -2944,10 +2959,68 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, GV, CS, p_atm) state%Hml(i,j) = depth(i) enddo enddo ! end of j loop - endif ! end BULKMIXEDLAYER - state%u => u(:,:,1) - state%v => v(:,:,1) + ! Determine the mean velocities in the uppermost depth_ml fluid. +!$OMP parallel do default(none) shared(is,ie,js,je,nz,CS,state,h,depth_ml,GV) private(depth,dh) + do j=js,je + + do i=is,ie + depth(i) = 0.0 + state%v(i,j) = 0.0 + hv(i,j,:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + enddo + + do k=1,nz ; do i=is,ie + if (depth(i) + hv(i,j,k)*GV%H_to_m < depth_ml) then + dh = hv(i,j,k)*GV%H_to_m + elseif (depth(i) < depth_ml) then + dh = depth_ml - depth(i) + else + dh = 0.0 + endif + state%v(i,j) = state%v(i,j) + dh * v(i,j,k) + depth(i) = depth(i) + dh + enddo ; enddo + ! Calculate the average properties of the mixed layer depth. + do i=is,ie + if (depth(i) < GV%H_subroundoff*GV%H_to_m) & + depth(i) = GV%H_subroundoff*GV%H_to_m + state%v(i,j) = state%v(i,j) / depth(i) + enddo + enddo ! end of j loop + +!$OMP parallel do default(none) shared(is,ie,js,je,nz,CS,state,h,depth_ml,GV) private(depth,dh) + do j=js,je + + do i=is,ie + depth(i) = 0.0 + state%u(i,j) = 0.0 + hu(i,j,:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + enddo + + do k=1,nz ; do i=is,ie + if (depth(i) + hu(i,j,k)*GV%H_to_m < depth_ml) then + dh = hu(i,j,k)*GV%H_to_m + elseif (depth(i) < depth_ml) then + dh = depth_ml - depth(i) + else + dh = 0.0 + endif + state%u(i,j) = state%u(i,j) + dh * u(i,j,k) + depth(i) = depth(i) + dh + enddo ; enddo + ! Calculate the average properties of the mixed layer depth. + do i=is,ie + if (depth(i) < GV%H_subroundoff*GV%H_to_m) & + depth(i) = GV%H_subroundoff*GV%H_to_m + state%u(i,j) = state%u(i,j) / depth(i) + enddo + enddo ! end of j loop + endif ! end BULKMIXEDLAYER + + !!!! + !state%u => u(:,:,1) + !state%v => v(:,:,1) state%frazil => CS%tv%frazil state%TempxPmE => CS%tv%TempxPmE state%internal_heat => CS%tv%internal_heat From f7d85a9c9e757d9282ff450c9cf5bc2d74a2397d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 31 Aug 2016 09:22:05 -0400 Subject: [PATCH 29/46] Added the option to average surface properties over a fixed distance and when using BULKMIXLAYER When BULKMIXLAYER is used, the logical USE_HMIX_SFC_PROP = True (default is False) controls whether the surface properties in the boundary layer are averged over a fixed distance (given by HMIX_SFC_PROP, default is 1 m) or values at k=1, i.e., state%SST => CS%tv%T(:,:,1), should be used. I've added this capability to check the sensitivity of melting rate under an ice shelf (using the ISOMIP setup) to HMIX_SFC_PROP. --- src/core/MOM.F90 | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 1542bc6ad5..05c908bad1 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -217,6 +217,9 @@ module MOM real :: Hmix !< diagnostic mixed layer thickness (meter) when !! bulk mixed layer is not used. + logical :: usehmix !< If true, surface properties (e.g., SST) are averaged + !! over a distance given by HMIX_SFC_PROP. This parameter + !! is valid just when BULKMIXLAYER is defined. real :: missing=-1.0e34 !< missing data value for masked fields ! Flags needed to reach between start and finish phases of initialization @@ -1533,6 +1536,23 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) "case is the largest integer multiple of the coupling \n"//& "timestep that is less than or equal to DT_THERM.", default=.false.) + if (CS%bulkmixedlayer) then + call get_param(param_file, "MOM", "USE_HMIX_SFC_PROP", CS%usehmix, & + "Applied only if BULKMIXLAYER is used. If TRUE, surface properties \n"//& + "like SST, SSS, density and velocities are averaged over \n"//& + "a distance given by HMIX_SFC_PROP. These averaged fields are \n"//& + "passed to the surface state structure, which is then used \n"//& + "to, for example, compute melting under an ice shelf.", default=.false.) + if(CS%usehmix) then + call get_param(param_file, "MOM", "HMIX_SFC_PROP", CS%Hmix, & + "If USE_HMIX_SFC_PROP is true, HMIX_SFC_PROP is the depth \n"//& + "over which to average to find surface properties like \n"//& + "SST, SSS, density are surface velocities.", & + units="m", default=1.0) + endif + + endif + if (.not.CS%bulkmixedlayer) then call get_param(param_file, "MOM", "HMIX_SFC_PROP", CS%Hmix, & "If BULKMIXEDLAYER is false, HMIX_SFC_PROP is the depth \n"//& @@ -2890,7 +2910,7 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, GV, CS, p_atm) enddo ; enddo endif ; endif - if (CS%bulkmixedlayer) then + if (CS%bulkmixedlayer .and. .not. CS%usehmix) then state%SST => CS%tv%T(:,:,1) state%SSS => CS%tv%S(:,:,1) state%u => u(:,:,1) From a2a25adb52f6bf8d8ab78a0043d7b49cf690f98e Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 31 Aug 2016 09:54:41 -0400 Subject: [PATCH 30/46] Added the option to compute S at the ice-ocean interface using a quadratic eq. By seting find_salt_root = True, Sbdry is computed using a quadratic equation (derived from the isomip 3 eq. formulation). If find_salt_root = False, the previous interactive method is used. I also fixed a bug in the ustar_shelf calculation. --- src/ice_shelf/MOM_ice_shelf.F90 | 123 ++++++++++++++++++++++---------- 1 file changed, 86 insertions(+), 37 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 9376cf8e60..78fed9fea8 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -336,8 +336,11 @@ module MOM_ice_shelf logical :: threeeq ! If true, the 3 equation consistency equations are ! used to calculate the flux at the ocean-ice interface. logical :: insulator ! If true, ice shelf is a perfect insulator - logical :: const_gamma ! If true, gamma_T is specified by the user. + logical :: const_gamma ! If true, gamma_T is specified by the user. + logical :: find_salt_root ! If true, if true find Sbdry using a quadratic eq. real :: cutoff_depth ! depth above which melt is set to zero (>= 0). + real :: lambda1, lambda2, lambda3 ! liquidus coeffs. Needed if + ! find_salt_root = true integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, & id_tfreeze = -1, id_tfl_shelf = -1, & id_thermal_driving = -1, id_haline_driving = -1, & @@ -449,8 +452,9 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) real :: PR, SC ! The Prandtl number and Schmidt number, nondim. ! 3 equation formulation variables - real :: Sbdry ! Salinities in the ocean at the interface with the - real :: Sbdry_it ! the ice shelf, in PSU. + real :: Sbdry ! Salinities in the ocean at the interface with the + real :: Sbdry_it ! the ice shelf, in PSU. + real :: Sbdry1, Sbdry2, S_a, S_b, S_c ! use to find salt roots real :: dS_it ! The interface salinity change during an iteration, in PSU. real :: hBL_neut ! The neutral boundary layer thickness, in m. real :: hBL_neut_h_molec ! The ratio of the neutral boundary layer thickness @@ -534,12 +538,12 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) ! Iteratively determine a self-consistent set of fluxes, with the ocean ! salinity just below the ice-shelf as the variable that is being ! iterated for. -! ### SHOULD I SET USTAR_SHELF YET? + ! ### SHOULD I SET USTAR_SHELF YET? u_at_h = state%u(i,j) v_at_h = state%v(i,j) - fluxes%ustar_shelf(i,j)= sqrt(CS%cdrag)*sqrt((u_at_h**2.0 + v_at_h**2.0)**2 +& + fluxes%ustar_shelf(i,j)= sqrt(CS%cdrag)*sqrt(sqrt(u_at_h**2.0 + v_at_h**2.0)**2 +& CS%utide(i,j)**2) ustar_h = MAX(CS%ustar_bg, fluxes%ustar_shelf(i,j)) @@ -559,14 +563,35 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) else ; hBL_neut = (VK*ustar_h) / absf ; endif hBL_neut_h_molec = ZETA_N * ((hBL_neut * ustar_h) / (5.0 * CS%Kv_molec)) - ! Guess sss as the iteration starting point for the boundary salinity. - Sbdry = state%sss(i,j) ; Sb_max_set = .false. ; Sb_min_set = .false. - ! Determine the mixed layer buoyancy flux, wB_flux. dB_dS = (CS%g_Earth / Rhoml(i)) * dR0_dS(i) dB_dT = (CS%g_Earth / Rhoml(i)) * dR0_dT(i) ln_neut = 0.0 ; if (hBL_neut_h_molec > 1.0) ln_neut = log(hBL_neut_h_molec) + if (CS%find_salt_root) then + ! read liquidus parameters + + S_a = CS%lambda1 * CS%Gamma_T_3EQ * CS%Cp +! S_b = -CS%Gamma_T_3EQ*(CS%lambda2-CS%lambda3*p_int(i)-state%sst(i,j)) & +! -LF*CS%Gamma_T_3EQ/35.0 + + S_b = CS%Gamma_T_3EQ*CS%Cp*(CS%lambda2+CS%lambda3*p_int(i)-state%sst(i,j)) & + -LF*CS%Gamma_T_3EQ/35.0 + S_c = LF*(CS%Gamma_T_3EQ/35.0)*state%sss(i,j) + !write(*,*)'a,b,c',S_a,S_b,S_c + + Sbdry1 = (-S_b + SQRT(S_b*S_b-4*S_a*S_c))/(2*S_a) + Sbdry2 = (-S_b - SQRT(S_b*S_b-4*S_a*S_c))/(2*S_a) + Sbdry = MAX(Sbdry1, Sbdry2) + !write(*,*)'#### Sbdry ####',Sbdry, Sbdry1, Sbdry2 + if (Sbdry < 0.) call MOM_error(FATAL, & + "shelf_calc_flux: Negative salinity (Sbdry).") + + else + ! Guess sss as the iteration starting point for the boundary salinity. + Sbdry = state%sss(i,j) ; Sb_max_set = .false. ; Sb_min_set = .false. + endif !find_salt_root + do it1 = 1,20 ! Determine the potential temperature at the ice-ocean interface. call calculate_TFreeze(Sbdry, p_int(i), CS%tfreeze(i,j), CS%eqn_of_state) @@ -644,7 +669,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) CS%exch_vel_t(i,j) = ustar_h * I_Gam_T CS%exch_vel_s(i,j) = ustar_h * I_Gam_S - !Calculate the heat flux inside the ice shelf. + !Calculate the heat flux inside the ice shelf. !vertical adv/diff as in H+J 1999, eqns (26) & approx from (31). ! Q_ice = rho_ice * CS%CP_Ice * K_ice * dT/dz (at interface) @@ -673,33 +698,39 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) endif endif - !other options: dTi/dz linear through shelf - ! dTi_dz = (CS%Temp_Ice - CS%tfreeze(i,j))/G%draft(i,j) - ! CS%tflux_shelf(i,j) = - Rho_Ice * CS%CP_Ice * KTI * dTi_dz - - mass_exch = CS%exch_vel_s(i,j) * CS%Rho0 - Sbdry_it = (state%sss(i,j) * mass_exch + CS%Salin_ice * CS%lprec(i,j)) / & - (mass_exch + CS%lprec(i,j)) - dS_it = Sbdry_it - Sbdry - if (abs(dS_it) < 1e-4*(0.5*(state%sss(i,j) + Sbdry + 1.e-10))) exit - - if (dS_it < 0.0) then ! Sbdry is now the upper bound. - if (Sb_max_set .and. (Sbdry > Sb_max)) call MOM_error(FATAL, & - "shelf_calc_flux: Irregular iteration for Sbdry (max).") - Sb_max = Sbdry ; dS_max = dS_it ; Sb_max_set = .true. - else ! Sbdry is now the lower bound. - if (Sb_min_set .and. (Sbdry < Sb_min)) call MOM_error(FATAL, & - "shelf_calc_flux: Irregular iteration for Sbdry (min).") - Sb_min = Sbdry ; dS_min = dS_it ; Sb_min_set = .true. + !other options: dTi/dz linear through shelf + ! dTi_dz = (CS%Temp_Ice - CS%tfreeze(i,j))/G%draft(i,j) + ! CS%tflux_shelf(i,j) = - Rho_Ice * CS%CP_Ice * KTI * dTi_dz + + + if (CS%find_salt_root) then + exit ! no need to do interaction, so exit loop + else + + mass_exch = CS%exch_vel_s(i,j) * CS%Rho0 + Sbdry_it = (state%sss(i,j) * mass_exch + CS%Salin_ice * CS%lprec(i,j)) / & + (mass_exch + CS%lprec(i,j)) + dS_it = Sbdry_it - Sbdry + if (abs(dS_it) < 1e-4*(0.5*(state%sss(i,j) + Sbdry + 1.e-10))) exit + + if (dS_it < 0.0) then ! Sbdry is now the upper bound. + if (Sb_max_set .and. (Sbdry > Sb_max)) call MOM_error(FATAL, & + "shelf_calc_flux: Irregular iteration for Sbdry (max).") + Sb_max = Sbdry ; dS_max = dS_it ; Sb_max_set = .true. + else ! Sbdry is now the lower bound. + if (Sb_min_set .and. (Sbdry < Sb_min)) call MOM_error(FATAL, & + "shelf_calc_flux: Irregular iteration for Sbdry (min).") + Sb_min = Sbdry ; dS_min = dS_it ; Sb_min_set = .true. + endif + if (Sb_min_set .and. Sb_max_set) then + ! Use the false position method for the next iteration. + Sbdry = Sb_min + (Sb_max-Sb_min) * (dS_min / (dS_min - dS_max)) + else + Sbdry = Sbdry_it + endif + + Sbdry = Sbdry_it endif - if (Sb_min_set .and. Sb_max_set) then - ! Use the false position method for the next iteration. - Sbdry = Sb_min + (Sb_max-Sb_min) * (dS_min / (dS_min - dS_max)) - else - Sbdry = Sbdry_it - endif - - Sbdry = Sbdry_it enddo !it1 ! Check for non-convergence and/or non-boundedness? @@ -1240,6 +1271,25 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti if (CS%const_gamma) call get_param(param_file, mod, "SHELF_3EQ_GAMMA_T", CS%Gamma_T_3EQ, & "Nondimensional heat-transfer coefficient.",default=2.2E-2, & units="nondim.", fail_if_missing=.true.) + if (CS%threeeq) & + call get_param(param_file, mod, "SHELF_S_ROOT", CS%find_salt_root, & + "If SHELF_S_ROOT = True, salinity at the ice/ocean interface (Sbdry) \n "//& + "is computed from a quadratic equation. Otherwise, the previous \n"//& + "interactive method to estimate Sbdry is used.", default=.false.) + if (CS%find_salt_root) then ! read liquidus coeffs. + call get_param(param_file, mod, "TFREEZE_S0_P0",CS%lambda1, & + "this is the freezing potential temperature at \n"//& + "S=0, P=0.", units="deg C", default=0.0, do_not_log=.true.) + call get_param(param_file, mod, "DTFREEZE_DS",CS%lambda1, & + "this is the derivative of the freezing potential \n"//& + "temperature with salinity.", & + units="deg C PSU-1", default=-0.054, do_not_log=.true.) + call get_param(param_file, mod, "DTFREEZE_DP",CS%lambda3, & + "this is the derivative of the freezing potential \n"//& + "temperature with pressure.", & + units="deg C Pa-1", default=0.0, do_not_log=.true.) + + endif if (.not.CS%threeeq) & call get_param(param_file, mod, "SHELF_2EQ_GAMMA_T", CS%gamma_t, & @@ -1424,7 +1474,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti allocate( CS%tflux_shelf(isd:ied,jsd:jed) ) ; CS%tflux_shelf(:,:) = 0.0 allocate( CS%tfreeze(isd:ied,jsd:jed) ) ; CS%tfreeze(:,:) = 0.0 -! allocate( CS%Sbdry(isd:ied,jsd:jed) ) ; CS%Sbdry(:,:) = 0.0 allocate( CS%exch_vel_s(isd:ied,jsd:jed) ) ; CS%exch_vel_s(:,:) = 0.0 allocate( CS%exch_vel_t(isd:ied,jsd:jed) ) ; CS%exch_vel_t(:,:) = 0.0 @@ -5626,7 +5675,7 @@ subroutine ice_shelf_end(CS) deallocate(CS%t_flux) ; deallocate(CS%lprec) deallocate(CS%salt_flux) - deallocate(CS%tflux_shelf) ; deallocate(CS%tfreeze);! deallocate(CS%Sbdry); + deallocate(CS%tflux_shelf) ; deallocate(CS%tfreeze); deallocate(CS%exch_vel_t) ; deallocate(CS%exch_vel_s) deallocate(CS%h_shelf) ; deallocate(CS%hmask) From 7c841215cd323181ddf210d18cb8158b36a4f76f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 31 Aug 2016 10:10:59 -0400 Subject: [PATCH 31/46] Flip the grid (x is longitude and y is latitude). --- src/user/ISOMIP_initialization.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 1bb1fa4b85..81f031c7a2 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -108,7 +108,7 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) if (is_2D) then do j=js,je ; do i=is,ie ! 2D setup - xtil = G%geoLatT(i,j)*1.0e3/xbar + xtil = G%geoLonT(i,j)*1.0e3/xbar !xtil = 450*1.0e3/xbar bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 !by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & @@ -127,18 +127,18 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) do j=js,je ; do i=is,ie ! 3D setup ! #### TEST ####### - !if (G%geoLatT(i,j)<500.) then + !if (G%geoLonT(i,j)<500.) then ! xtil = 500.*1.0e3/xbar !else - ! xtil = G%geoLatT(i,j)*1.0e3/xbar + ! xtil = G%geoLonT(i,j)*1.0e3/xbar !endif ! ##### TEST ##### - xtil = G%geoLatT(i,j)*1.0e3/xbar + xtil = G%geoLonT(i,j)*1.0e3/xbar bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 - by = (dc/(1.+exp(-2.*(G%geoLonT(i,j)*1.0e3- ly/2. - wc)/fc))) + & - (dc/(1.+exp(2.*(G%geoLonT(i,j)*1.0e3- ly/2. + wc)/fc))) + by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & + (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) D(i,j) = -max(bx+by,-bmax) if (D(i,j) > max_depth) D(i,j) = max_depth From 9cec7dff9d7e44b9da2874905856c2f022f4e852 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 31 Aug 2016 10:18:10 -0400 Subject: [PATCH 32/46] Replace geoLatT to geoLonT in the sponge subroutine. --- src/user/ISOMIP_initialization.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 81f031c7a2..1588038df1 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -486,10 +486,10 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) ! and mask2dT is 1. do i=is,ie; do j=js,je - if (G%geoLatT(i,j) >= 790.0 .AND. G%geoLatT(i,j) <= 800.0) then + if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then ! 1 / day - dummy1=(G%geoLatT(i,j)-790.0)/(800.0-790.0) + dummy1=(G%geoLonT(i,j)-790.0)/(800.0-790.0) damp = 1.0/TNUDG * max(0.0,dummy1) else ; damp=0.0 From 7939f7eaac25df59cce52c7e7572c081936737ba Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 1 Sep 2016 10:15:00 -0400 Subject: [PATCH 33/46] Added some print statments for debugging --- src/core/MOM.F90 | 2 -- src/ice_shelf/MOM_ice_shelf.F90 | 9 +++++++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 05c908bad1..7e9234630f 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2900,7 +2900,6 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, GV, CS, p_atm) iscB = G%iscB ; iecB = G%iecB; jscB = G%jscB ; jecB = G%jecB isdB = G%isdB ; iedB = G%iedB; jsdB = G%jsdB ; jedB = G%jedB - !write(*,*)'iscB,iecB,jscB,jecB', iscB,iecB,jscB,jecB state%sea_lev => ssh if (present(p_atm)) then ; if (ASSOCIATED(p_atm)) then @@ -3042,7 +3041,6 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, GV, CS, p_atm) enddo ! end of j loop endif ! end BULKMIXEDLAYER - !!!! !state%u => u(:,:,1) !state%v => v(:,:,1) state%frazil => CS%tv%frazil diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 78fed9fea8..05835d2aca 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -584,9 +584,14 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) Sbdry2 = (-S_b - SQRT(S_b*S_b-4*S_a*S_c))/(2*S_a) Sbdry = MAX(Sbdry1, Sbdry2) !write(*,*)'#### Sbdry ####',Sbdry, Sbdry1, Sbdry2 - if (Sbdry < 0.) call MOM_error(FATAL, & + if (Sbdry < 0.) then + ! this is for debuging! + write(*,*)'state%sss(i,j)',state%sss(i,j) + write(*,*)'S_a, S_b, S_c',S_a, S_b, S_c + write(*,*)'I,J,Sbdry1,Sbdry2',i,j,Sbdry1,Sbdry2 + call MOM_error(FATAL, & "shelf_calc_flux: Negative salinity (Sbdry).") - + endif else ! Guess sss as the iteration starting point for the boundary salinity. Sbdry = state%sss(i,j) ; Sb_max_set = .false. ; Sb_min_set = .false. From df968c9d8f1f9499ed0b886556d004dcd96cb233 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 2 Sep 2016 15:18:49 -0400 Subject: [PATCH 34/46] Added the capability to remap to a z* coordinate underneath an ice shelf --- src/core/MOM.F90 | 2 +- src/diagnostics/MOM_diag_to_Z.F90 | 295 ++++++++++++++++++------------ 2 files changed, 177 insertions(+), 120 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 7e9234630f..a9dc54a668 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1161,7 +1161,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) if (Time_local + set_time(int(0.5*dt_therm)) > CS%Z_diag_time) then call enable_averaging(real(time_type_to_real(CS%Z_diag_interval)), & CS%Z_diag_time, CS%diag) - call calculate_Z_diag_fields(u, v, h, CS%dt_trans, & + call calculate_Z_diag_fields(u, v, h, ssh,CS%dt_trans, & G, GV, CS%diag_to_Z_CSp) CS%Z_diag_time = CS%Z_diag_time + CS%Z_diag_interval call disable_averaging(CS%diag) diff --git a/src/diagnostics/MOM_diag_to_Z.F90 b/src/diagnostics/MOM_diag_to_Z.F90 index 4c019bd9d8..9c0a1744bb 100644 --- a/src/diagnostics/MOM_diag_to_Z.F90 +++ b/src/diagnostics/MOM_diag_to_Z.F90 @@ -42,7 +42,7 @@ module MOM_diag_to_Z !* The boundaries always run through q grid points (x). * !* * !********+*********+*********+*********+*********+*********+*********+** - +use MOM_domains, only : pass_var use MOM_coms, only : reproducing_sum use MOM_diag_mediator, only : post_data, post_data_1d_k, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl, time_type, diag_axis_init @@ -162,27 +162,30 @@ function global_z_mean(var,G,CS,tracer) end function global_z_mean -subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) +subroutine calculate_Z_diag_fields(u, v, h, ssh_in, dt, G, GV, CS) type(ocean_grid_type), intent(inout) :: G type(verticalGrid_type), intent(in) :: GV real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_in real, intent(in) :: dt type(diag_to_Z_CS), pointer :: CS ! This subroutine maps tracers and velocities into depth space for diagnostics. ! Arguments: -! (in) u - zonal velocity component (m/s) -! (in) v - meridional velocity component (m/s) -! (in) h - layer thickness (meter or kg/m2) -! (in) dt - time difference in sec since last call to this subroutine -! (in) G - ocean grid structure -! (in) GV - The ocean's vertical grid structure. -! (in) CS - control structure returned by previous call to diagnostics_init +! (in) u - zonal velocity component (m/s) +! (in) v - meridional velocity component (m/s) +! (in) h - layer thickness (meter or kg/m2) +! (in) ssh_in - sea surface height (meter or kg/m2) +! (in) dt - time difference in sec since last call to this subroutine +! (in) G - ocean grid structure +! (in) GV - The ocean's vertical grid structure. +! (in) CS - control structure returned by previous call to diagnostics_init ! Note the deliberately reversed axes in h_f, u_f, v_f, and tr_f. + real :: ssh(SZI_(G),SZJ_(G)) ! copy of ssh_in (meter or kg/m2) real :: e(SZK_(G)+2) ! z-star interface heights (meter or kg/m2) real :: h_f(SZK_(G)+1,SZI_(G)) ! thicknesses of massive layers (meter or kg/m2) real :: u_f(SZK_(G)+1,SZIB_(G))! zonal velocity component in any massive layer @@ -191,7 +194,8 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) real :: tr_f(SZK_(G),max(CS%num_tr_used,1),SZI_(G)) ! tracer concentration in massive layers integer :: nk_valid(SZIB_(G)) ! number of massive layers in a column - real :: D_pt(SZIB_(G)) ! bottom depth (meter or kg/m2) + real :: D_pt(SZIB_(G)) ! bottom depth (meter or kg/m2) + real :: shelf_depth(SZIB_(G)) ! ice shelf depth (meter or kg/m2) real :: htot ! summed layer thicknesses (meter or kg/m2) real :: dilate ! proportion by which to dilate every layer real :: wt(SZK_(G)+1) ! fractional weight for each layer in the @@ -212,11 +216,15 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) integer :: k_top, k_bot, k_bot_prev integer :: i, j, k, k2, kz, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nk, m, nkml - + integer :: IsgB, IegB, JsgB, JegB is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nk = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + IsgB = G%IsgB ; IegB = G%IegB ; JsgB = G%JsgB ; JegB = G%JegB nkml = max(GV%nkml, 1) Angstrom = GV%Angstrom + ssh(:,:) = ssh_in + ! Update the halos + call pass_var(ssh, G%Domain) if (.not.associated(CS)) call MOM_error(FATAL, & "diagnostic_fields_zstar: Module must be initialized before it is used.") @@ -231,10 +239,21 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) CS%u_z(I,j,kz) = CS%missing_vel enddo ; enddo ; enddo + do j=js,je ! Remove all massless layers. do I=Isq,Ieq - nk_valid(I) = 0 ; D_pt(I) = 0.5*(G%bathyT(i+1,j)+G%bathyT(i,j)) + nk_valid(I) = 0 + D_pt(I) = 0.5*(G%bathyT(i+1,j)+G%bathyT(i,j)) + ! GM: The following is workaround to make this subroutine works under an ice shelf + ! Zero ssh in the open ocean and get the depth of ice shelf (hence the abs below) + ! I am not sure if -0.1 is the best choice + shelf_depth(I) = 0.5*(ssh(i+1,j)+ssh(i,j)) + if (shelf_depth(I) > -0.1) then ! open ocean + shelf_depth(I) = 0.0 + else ! ice shelf + shelf_depth(I) = abs(shelf_depth(I)) + endif enddo do k=1,nk ; do I=Isq,Ieq if ((G%mask2dCu(I,j) > 0.5) .and. (h(i,j,k)+h(i+1,j,k) > 4.0*Angstrom)) then @@ -247,13 +266,22 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) ! no-slip BBC in the output, if anything but piecewise constant is used. nk_valid(I) = nk_valid(I) + 1 ; k2 = nk_valid(I) h_f(k2,I) = Angstrom ; u_f(k2,I) = 0.0 + ! GM: D_pt is always slightly larger (by 1E-6 or so) than shelf_depth, so + ! I consider that the ice shelf is grounded when + ! shelf_depth(I) + 1.0E-3 > D_pt(i) + if (shelf_depth(I) + 1.0E-3 > D_pt(i)) nk_valid(I)=0 endif ; enddo + do I=Isq,Ieq ; if (nk_valid(I) > 0) then ! Calculate the z* interface heights for tracers. htot = 0.0 ; do k=1,nk_valid(i) ; htot = htot + h_f(k,i) ; enddo - dilate = 0.0 ; if (htot*GV%H_to_m > 2.0*Angstrom) dilate = (D_pt(i) - 0.0) / htot - + dilate = 0.0 + if (htot*GV%H_to_m > 2.0*Angstrom) then + dilate = MAX((D_pt(i) - shelf_depth(i)),Angstrom)/htot + !dilate = (D_pt(i) - 0.0)/htot + !write(*,*)MAX((D_pt(i) - shelf_depth(i)),Angstrom)/htot + endif e(nk_valid(i)+1) = -D_pt(i) do k=nk_valid(i),1,-1 ; e(K) = e(K+1) + h_f(k,i)*dilate ; enddo @@ -264,39 +292,44 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) k_bot, k_top, k_bot, wt, z1, z2) if (k_top>nk_valid(I)) exit - if (linear_velocity_profiles) then - k = k_top - if (k /= k_bot_prev) then - ! Calculate the intra-cell profile. - slope = 0.0 ! ; curv = 0.0 - if ((k < nk_valid(I)) .and. (k > nkml)) call & - find_limited_slope(u_f(:,I), e, slope, k) - endif - ! This is the piecewise linear form. - CS%u_z(I,j,kz) = wt(k) * (u_f(k,I) + 0.5*slope*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - do k=k_top+1,k_bot-1 - CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k)*u_f(k,I) - enddo - if (k_bot > k_top) then ; k = k_bot - ! Calculate the intra-cell profile. - slope = 0.0 ! ; curv = 0.0 - if ((k < nk_valid(I)) .and. (k > nkml)) call & - find_limited_slope(u_f(:,I), e, slope, k) - ! This is the piecewise linear form. - CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k) * & - (u_f(k,I) + 0.5*slope*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - endif - k_bot_prev = k_bot - else ! Use piecewise constant profiles. - CS%u_z(I,j,kz) = wt(k_top)*u_f(k_top,I) - do k=k_top+1,k_bot - CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k)*u_f(k,I) - enddo - endif ! linear profiles + !GM if top range that is being map is below the shelf, interpolate + ! otherwise keep missing_vel + if (CS%Z_int(kz)<=-shelf_depth(I)) then + + if (linear_velocity_profiles) then + k = k_top + if (k /= k_bot_prev) then + ! Calculate the intra-cell profile. + slope = 0.0 ! ; curv = 0.0 + if ((k < nk_valid(I)) .and. (k > nkml)) call & + find_limited_slope(u_f(:,I), e, slope, k) + endif + ! This is the piecewise linear form. + CS%u_z(I,j,kz) = wt(k) * (u_f(k,I) + 0.5*slope*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + do k=k_top+1,k_bot-1 + CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k)*u_f(k,I) + enddo + if (k_bot > k_top) then ; k = k_bot + ! Calculate the intra-cell profile. + slope = 0.0 ! ; curv = 0.0 + if ((k < nk_valid(I)) .and. (k > nkml)) call & + find_limited_slope(u_f(:,I), e, slope, k) + ! This is the piecewise linear form. + CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k) * & + (u_f(k,I) + 0.5*slope*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + endif + k_bot_prev = k_bot + else ! Use piecewise constant profiles. + CS%u_z(I,j,kz) = wt(k_top)*u_f(k_top,I) + do k=k_top+1,k_bot + CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k)*u_f(k,I) + enddo + endif ! linear profiles + endif ! below shelf enddo ! kz-loop endif ; enddo ! I-loop and mask enddo ! j-loop @@ -314,6 +347,12 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) ! Remove all massless layers. do i=is,ie nk_valid(i) = 0 ; D_pt(i) = 0.5*(G%bathyT(i,j)+G%bathyT(i,j+1)) + shelf_depth(I) = 0.5*(ssh(i,j)+ssh(i,j+1)) + if (shelf_depth(I) > -0.1) then ! open ocean + shelf_depth(I) = 0.0 + else ! ice shelf + shelf_depth(I) = abs(shelf_depth(I)) + endif enddo do k=1,nk ; do i=is,ie if ((G%mask2dCv(i,j) > 0.5) .and. (h(i,j,k)+h(i,j+1,k) > 4.0*Angstrom)) then @@ -326,13 +365,16 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) ! no-slip BBC in the output, if anything but piecewise constant is used. nk_valid(i) = nk_valid(i) + 1 ; k2 = nk_valid(i) h_f(k2,i) = Angstrom ; v_f(k2,i) = 0.0 + if (shelf_depth(I) + 1.0E-3 > D_pt(i)) nk_valid(I)=0 endif ; enddo do i=is,ie ; if (nk_valid(i) > 0) then ! Calculate the z* interface heights for tracers. htot = 0.0 ; do k=1,nk_valid(i) ; htot = htot + h_f(k,i) ; enddo - dilate = 0.0 ; if (htot > 2.0*Angstrom) dilate = (D_pt(i) - 0.0) / htot - + dilate = 0.0 + if (htot > 2.0*Angstrom) then + dilate = MAX((D_pt(i) - shelf_depth(i)),Angstrom)/htot + endif e(nk_valid(i)+1) = -D_pt(i) do k=nk_valid(i),1,-1 ; e(K) = e(K+1) + h_f(k,i)*dilate ; enddo @@ -342,41 +384,44 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) call find_overlap(e, CS%Z_int(kz), CS%Z_int(kz+1), nk_valid(i), & k_bot, k_top, k_bot, wt, z1, z2) if (k_top>nk_valid(i)) exit - - if (linear_velocity_profiles) then - k = k_top - if (k /= k_bot_prev) then - ! Calculate the intra-cell profile. - slope = 0.0 ! ; curv = 0.0 - if ((k < nk_valid(i)) .and. (k > nkml)) call & - find_limited_slope(v_f(:,i), e, slope, k) - endif - ! This is the piecewise linear form. - CS%v_z(i,J,kz) = wt(k) * (v_f(k,i) + 0.5*slope*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - do k=k_top+1,k_bot-1 - CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k)*v_f(k,i) - enddo - if (k_bot > k_top) then ; k = k_bot - ! Calculate the intra-cell profile. - slope = 0.0 ! ; curv = 0.0 - if ((k < nk_valid(i)) .and. (k > nkml)) call & - find_limited_slope(v_f(:,i), e, slope, k) - ! This is the piecewise linear form. - CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k) * & - (v_f(k,i) + 0.5*slope*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - endif - k_bot_prev = k_bot - else ! Use piecewise constant profiles. - CS%v_z(i,J,kz) = wt(k_top)*v_f(k_top,i) - do k=k_top+1,k_bot - CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k)*v_f(k,i) - enddo - endif ! linear profiles - enddo ! kz-loop + !GM if top range that is being map is below the shelf, interpolate + ! otherwise keep missing_vel + if (CS%Z_int(kz)<=-shelf_depth(I)) then + if (linear_velocity_profiles) then + k = k_top + if (k /= k_bot_prev) then + ! Calculate the intra-cell profile. + slope = 0.0 ! ; curv = 0.0 + if ((k < nk_valid(i)) .and. (k > nkml)) call & + find_limited_slope(v_f(:,i), e, slope, k) + endif + ! This is the piecewise linear form. + CS%v_z(i,J,kz) = wt(k) * (v_f(k,i) + 0.5*slope*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + do k=k_top+1,k_bot-1 + CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k)*v_f(k,i) + enddo + if (k_bot > k_top) then ; k = k_bot + ! Calculate the intra-cell profile. + slope = 0.0 ! ; curv = 0.0 + if ((k < nk_valid(i)) .and. (k > nkml)) call & + find_limited_slope(v_f(:,i), e, slope, k) + ! This is the piecewise linear form. + CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k) * & + (v_f(k,i) + 0.5*slope*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + endif + k_bot_prev = k_bot + else ! Use piecewise constant profiles. + CS%v_z(i,J,kz) = wt(k_top)*v_f(k_top,i) + do k=k_top+1,k_bot + CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k)*v_f(k,i) + enddo + endif ! linear profiles + endif ! below shelf + enddo ! kz-loop endif ; enddo ! i-loop and mask enddo ! J-loop @@ -392,11 +437,20 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) do j=js,je ! Remove all massless layers. - do i=is,ie ; nk_valid(i) = 0 ; D_pt(i) = G%bathyT(i,j) ; enddo + do i=is,ie + nk_valid(i) = 0 ; D_pt(i) = G%bathyT(i,j) + shelf_depth(I) = ssh(i,j) + if (shelf_depth(I) > -0.1) then ! open ocean + shelf_depth(I) = 0.0 + else ! ice shelf + shelf_depth(I) = abs(shelf_depth(I)) + endif + enddo do k=1,nk ; do i=is,ie if ((G%mask2dT(i,j) > 0.5) .and. (h(i,j,k) > 2.0*Angstrom)) then nk_valid(i) = nk_valid(i) + 1 ; k2 = nk_valid(i) h_f(k2,i) = h(i,j,k) + if (shelf_depth(I) + 1.0E-3 > D_pt(i)) nk_valid(I)=0 do m=1,CS%num_tr_used ; tr_f(k2,m,i) = CS%tr_model(m)%p(i,j,k) ; enddo endif enddo ; enddo @@ -404,8 +458,10 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) do i=is,ie ; if (nk_valid(i) > 0) then ! Calculate the z* interface heights for tracers. htot = 0.0 ; do k=1,nk_valid(i) ; htot = htot + h_f(k,i) ; enddo - dilate = 0.0 ; if (htot > 2.0*Angstrom) dilate = (D_pt(i) - 0.0) / htot - + dilate = 0.0 + if (htot > 2.0*Angstrom) then + dilate = MAX((D_pt(i) - shelf_depth(i)),Angstrom)/htot + endif e(nk_valid(i)+1) = -D_pt(i) do k=nk_valid(i),1,-1 ; e(K) = e(K+1) + h_f(k,i)*dilate ; enddo @@ -415,38 +471,39 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) call find_overlap(e, CS%Z_int(kz), CS%Z_int(kz+1), nk_valid(i), & k_bot, k_top, k_bot, wt, z1, z2) if (k_top>nk_valid(i)) exit - - do m=1,CS%num_tr_used - k = k_top - if (k /= k_bot_prev) then - ! Calculate the intra-cell profile. - sl_tr(m) = 0.0 ! ; cur_tr(m) = 0.0 - if ((k < nk_valid(i)) .and. (k > nkml)) call & - find_limited_slope(tr_f(:,m,i), e, sl_tr(m), k) - endif - ! This is the piecewise linear form. - CS%tr_z(m)%p(i,j,kz) = wt(k) * & - (tr_f(k,m,i) + 0.5*sl_tr(m)*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*cur_tr(m)*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - do k=k_top+1,k_bot-1 - CS%tr_z(m)%p(i,j,kz) = CS%tr_z(m)%p(i,j,kz) + wt(k)*tr_f(k,m,i) - enddo - if (k_bot > k_top) then - k = k_bot - ! Calculate the intra-cell profile. - sl_tr(m) = 0.0 ! ; cur_tr(m) = 0.0 - if ((k < nk_valid(i)) .and. (k > nkml)) call & - find_limited_slope(tr_f(:,m,i), e, sl_tr(m), k) - ! This is the piecewise linear form. - CS%tr_z(m)%p(i,j,kz) = CS%tr_z(m)%p(i,j,kz) + wt(k) * & - (tr_f(k,m,i) + 0.5*sl_tr(m)*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*cur_tr(m)*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - endif - enddo - k_bot_prev = k_bot - enddo ! kz-loop + if (CS%Z_int(kz)<=-shelf_depth(I)) then + do m=1,CS%num_tr_used + k = k_top + if (k /= k_bot_prev) then + ! Calculate the intra-cell profile. + sl_tr(m) = 0.0 ! ; cur_tr(m) = 0.0 + if ((k < nk_valid(i)) .and. (k > nkml)) call & + find_limited_slope(tr_f(:,m,i), e, sl_tr(m), k) + endif + ! This is the piecewise linear form. + CS%tr_z(m)%p(i,j,kz) = wt(k) * & + (tr_f(k,m,i) + 0.5*sl_tr(m)*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*cur_tr(m)*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + do k=k_top+1,k_bot-1 + CS%tr_z(m)%p(i,j,kz) = CS%tr_z(m)%p(i,j,kz) + wt(k)*tr_f(k,m,i) + enddo + if (k_bot > k_top) then + k = k_bot + ! Calculate the intra-cell profile. + sl_tr(m) = 0.0 ! ; cur_tr(m) = 0.0 + if ((k < nk_valid(i)) .and. (k > nkml)) call & + find_limited_slope(tr_f(:,m,i), e, sl_tr(m), k) + ! This is the piecewise linear form. + CS%tr_z(m)%p(i,j,kz) = CS%tr_z(m)%p(i,j,kz) + wt(k) * & + (tr_f(k,m,i) + 0.5*sl_tr(m)*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*cur_tr(m)*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + endif + enddo + k_bot_prev = k_bot + endif ! below shelf + enddo ! kz-loop endif ; enddo ! i-loop and mask enddo ! j-loop From 028c52358be11d4c8e4740399f98a21bd51b61c4 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 5 Sep 2016 12:41:21 -0400 Subject: [PATCH 35/46] Normalize the amount of tracer based on max(melt) The concentration of passive tracer injected under the ice shelf is noe normalized by the maximum melting at that time. If melt > 0, the normalized tracer is then injected in the first two layers (ML). --- src/tracer/ISOMIP_tracer.F90 | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index 027e06389d..88bca36b73 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -33,7 +33,7 @@ module ISOMIP_tracer use MOM_variables, only : surface use MOM_open_boundary, only : ocean_OBC_type use MOM_verticalGrid, only : verticalGrid_type - +use MOM_coms, only : max_across_PEs use coupler_util, only : set_coupler_values, ind_csurf use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux @@ -333,19 +333,27 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G real :: mmax real :: b1(SZI_(G)) ! b1 and c1 are variables used by the real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver. + real :: melt(SZI_(G),SZJ_(G)) ! melt water (positive for melting + ! negative for freezing) integer :: i, j, k, is, ie, js, je, nz, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke if (.not.associated(CS)) return + melt(:,:) = fluxes%iceshelf_melt ! max. melt - mmax = MAXVAL(fluxes%iceshelf_melt) + mmax = MAXVAL(melt(is:ie,js:je)) + call max_across_PEs(mmax) + !write(*,*)'max melt', mmax ! dye melt water (m=1), dye = 1 if melt=max(melt) do m=1,NTR do j=js,je ; do i=is,ie - if (fluxes%iceshelf_melt(i,j) > 0.0) then - CS%tr(i,j,1,m) = fluxes%iceshelf_melt(i,j)/mmax + if (melt(i,j) > 0.0) then ! melting + !write(*,*)'i,j,melt,melt/mmax',i,j,melt(i,j),melt(i,j)/mmax + CS%tr(i,j,1:2,m) = melt(i,j)/mmax ! inject dye in the ML + else ! freezing + CS%tr(i,j,1:2,m) = 0.0 endif enddo ; enddo enddo From d33cc7ef93863addc3ca0b70b58fd67f8fc21ddc Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 9 Sep 2016 15:58:48 -0400 Subject: [PATCH 36/46] Fixed a bug in the ustar calculation. --- src/ice_shelf/MOM_ice_shelf.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 05835d2aca..c4d6cd83d2 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -543,8 +543,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) u_at_h = state%u(i,j) v_at_h = state%v(i,j) - fluxes%ustar_shelf(i,j)= sqrt(CS%cdrag)*sqrt(sqrt(u_at_h**2.0 + v_at_h**2.0)**2 +& - CS%utide(i,j)**2) + fluxes%ustar_shelf(i,j)= sqrt(CS%cdrag*(sqrt(u_at_h**2.0 + v_at_h**2.0)**2 +& + CS%utide(i,j)**2)) ustar_h = MAX(CS%ustar_bg, fluxes%ustar_shelf(i,j)) fluxes%ustar_shelf(i,j) = ustar_h From bbe45f5459a5e379f5352bb20df398763b96ac22 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 9 Sep 2016 23:28:32 -0400 Subject: [PATCH 37/46] Changed sponge code for layer mode I now read two different files, one with temp. and one with salt. values. This is work around to avoid having wrong T or S (depending on FIT_SALINITY) values near the surface within the sponge layer. To get salt values right in the sponge, FIT_SALINITY must be False. The oposite is true for temp. I could combine the *correct* temp and salt values in one file instead, but I prefer making this modification for now. --- src/user/ISOMIP_initialization.F90 | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 1588038df1..da1bbf7c22 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -437,8 +437,8 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) ! positive upward, in m. real :: min_depth, dummy1, z, delta_h real :: damp, rho_dummy, min_thickness, rho_tmp, xi0 - character(len=40) :: verticalCoordinate, filename, state_file, inputdir - character(len=40) :: temp_var, salt_var, eta_var + character(len=40) :: verticalCoordinate, filename, state_file_t, state_file_s + character(len=40) :: temp_var, salt_var, eta_var, inputdir character(len=40) :: mod = "ISOMIP_initialize_sponges" ! This subroutine's name. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz @@ -604,9 +604,17 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) ! 1) Read eta, salt and temp from IC file call get_param(PF, mod, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) - call get_param(PF, mod, "ISOMIP_SPONGE_FILE", state_file, & - "The name of the file with the state to damp toward.", & - fail_if_missing=.true.) + ! GM: get twp different files, one with temp and one with salt values + ! this is work around to avoid having wrong values near the surface + ! because of the FIT_SALINITY option. To get salt values right in the + ! sponge, FIT_SALINITY=False. The oposite is true for temp. One can + ! combined the *correct* temp and salt values in one file instead. + call get_param(PF, mod, "ISOMIP_SPONGE_FILE_TEMP", state_file_t, & + "The name of the file with temperatures and interfaces to \n"// & + " damp toward.", fail_if_missing=.true.) + call get_param(PF, mod, "ISOMIP_SPONGE_FILE_SALT", state_file_s, & + "The name of the file with salinities and interfaces to \n"// & + " damp toward.", fail_if_missing=.true.) call get_param(PF, mod, "SPONGE_PTEMP_VAR", temp_var, & "The name of the potential temperature variable in \n"//& @@ -618,14 +626,17 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) "The name of the interface height variable in \n"//& "SPONGE_STATE_FILE.", default="eta") - filename = trim(inputdir)//trim(state_file) - + !read temp and eta + filename = trim(inputdir)//trim(state_file_t) if (.not.file_exists(filename, G%Domain)) & call MOM_error(FATAL, " ISOMIP_initialize_sponges: Unable to open "//trim(filename)) - - !call read_data("",eta_var,eta(:,:,:), domain=G%Domain%mpp_domain) call read_data(filename,eta_var,eta(:,:,:), domain=G%Domain%mpp_domain) call read_data(filename,temp_var,T(:,:,:), domain=G%Domain%mpp_domain) + + ! read salt + filename = trim(inputdir)//trim(state_file_s) + if (.not.file_exists(filename, G%Domain)) & + call MOM_error(FATAL, " ISOMIP_initialize_sponges: Unable to open "//trim(filename)) call read_data(filename,salt_var,S(:,:,:), domain=G%Domain%mpp_domain) ! for debugging From 1114b2a813b8fde3c5ae42f0aab7b769ce0e610f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Sat, 10 Sep 2016 17:42:50 -0400 Subject: [PATCH 38/46] Going back to have just one file tha specifies the sponge layer properties in layer mode. --- src/user/ISOMIP_initialization.F90 | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index da1bbf7c22..670258c9ba 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -437,7 +437,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) ! positive upward, in m. real :: min_depth, dummy1, z, delta_h real :: damp, rho_dummy, min_thickness, rho_tmp, xi0 - character(len=40) :: verticalCoordinate, filename, state_file_t, state_file_s + character(len=40) :: verticalCoordinate, filename, state_file character(len=40) :: temp_var, salt_var, eta_var, inputdir character(len=40) :: mod = "ISOMIP_initialize_sponges" ! This subroutine's name. @@ -604,18 +604,14 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) ! 1) Read eta, salt and temp from IC file call get_param(PF, mod, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) - ! GM: get twp different files, one with temp and one with salt values + ! GM: get two different files, one with temp and one with salt values ! this is work around to avoid having wrong values near the surface ! because of the FIT_SALINITY option. To get salt values right in the ! sponge, FIT_SALINITY=False. The oposite is true for temp. One can ! combined the *correct* temp and salt values in one file instead. - call get_param(PF, mod, "ISOMIP_SPONGE_FILE_TEMP", state_file_t, & - "The name of the file with temperatures and interfaces to \n"// & + call get_param(PF, mod, "ISOMIP_SPONGE_FILE", state_file, & + "The name of the file with temps., salts. and interfaces to \n"// & " damp toward.", fail_if_missing=.true.) - call get_param(PF, mod, "ISOMIP_SPONGE_FILE_SALT", state_file_s, & - "The name of the file with salinities and interfaces to \n"// & - " damp toward.", fail_if_missing=.true.) - call get_param(PF, mod, "SPONGE_PTEMP_VAR", temp_var, & "The name of the potential temperature variable in \n"//& "SPONGE_STATE_FILE.", default="Temp") @@ -627,16 +623,11 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) "SPONGE_STATE_FILE.", default="eta") !read temp and eta - filename = trim(inputdir)//trim(state_file_t) + filename = trim(inputdir)//trim(state_file) if (.not.file_exists(filename, G%Domain)) & call MOM_error(FATAL, " ISOMIP_initialize_sponges: Unable to open "//trim(filename)) call read_data(filename,eta_var,eta(:,:,:), domain=G%Domain%mpp_domain) call read_data(filename,temp_var,T(:,:,:), domain=G%Domain%mpp_domain) - - ! read salt - filename = trim(inputdir)//trim(state_file_s) - if (.not.file_exists(filename, G%Domain)) & - call MOM_error(FATAL, " ISOMIP_initialize_sponges: Unable to open "//trim(filename)) call read_data(filename,salt_var,S(:,:,:), domain=G%Domain%mpp_domain) ! for debugging From 252098f17aec8932347376cabbd47597a7ea82b1 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 16 Sep 2016 14:43:56 -0400 Subject: [PATCH 39/46] Coded ustar in a cleaner way --- src/ice_shelf/MOM_ice_shelf.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index c4d6cd83d2..7cc1ecd724 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -543,7 +543,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) u_at_h = state%u(i,j) v_at_h = state%v(i,j) - fluxes%ustar_shelf(i,j)= sqrt(CS%cdrag*(sqrt(u_at_h**2.0 + v_at_h**2.0)**2 +& + fluxes%ustar_shelf(i,j)= sqrt(CS%cdrag*((u_at_h**2.0 + v_at_h**2.0) +& CS%utide(i,j)**2)) ustar_h = MAX(CS%ustar_bg, fluxes%ustar_shelf(i,j)) From 62f92a460986553d5c64d25301e70dc5150babee Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 20 Sep 2016 16:11:02 -0400 Subject: [PATCH 40/46] Added hML to the RESTART file and cleaned the ice shelf module hML is needed when using a restart file with the ice shelf module, so I added this field on MOM.F90 I added addional hchksum in MOM_ice_shelf.F9 to facilitate debudding. --- src/core/MOM.F90 | 6 ++++++ src/ice_shelf/MOM_ice_shelf.F90 | 37 +++++++++++++++++++++++---------- 2 files changed, 32 insertions(+), 11 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index a9dc54a668..67903aaeb5 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2861,6 +2861,12 @@ subroutine set_restart_fields(GV, param_file, CS) vd = var_desc("ave_ssh","meter","Time average sea surface height",'h','1') call register_restart_field(CS%ave_ssh, vd, .false., CS%restart_CSp) + ! GM, hML is needed when using the ice shelf module + if (associated(CS%tv%Hml)) then + vd = var_desc("hML","meter","Mixed layer thickness",'h','1') + call register_restart_field(CS%tv%Hml, vd, .false., CS%restart_CSp) + endif + end subroutine set_restart_fields diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 7cc1ecd724..6ffc8cde90 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -545,6 +545,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) fluxes%ustar_shelf(i,j)= sqrt(CS%cdrag*((u_at_h**2.0 + v_at_h**2.0) +& CS%utide(i,j)**2)) + + ustar_h = MAX(CS%ustar_bg, fluxes%ustar_shelf(i,j)) fluxes%ustar_shelf(i,j) = ustar_h @@ -784,7 +786,9 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) ALLOCATE ( mass_flux(G%ied,G%jed) ) mass_flux = (CS%lprec) * CS%area_shelf_h if (CS%DEBUG) then - call hchksum (CS%h_shelf, "melt rate", G%HI, haloshift=0) + call hchksum (fluxes%iceshelf_melt, "melt rate", G%HI, haloshift=0) + call hchksum (fluxes%ustar_shelf, "ustar_shelf calc", G%HI, haloshift=0) + call hchksum (state%Hml, "Hml", G%HI, haloshift=0) endif if (CS%shelf_mass_is_dynamic) then @@ -793,7 +797,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) call pass_var(CS%mass_shelf, G%domain) call cpu_clock_end(id_clock_pass) endif - + call add_shelf_flux(G, CS, state, fluxes) ! now the thermodynamic data is passed on... time to update the ice dynamic quantities @@ -923,6 +927,7 @@ subroutine add_shelf_flux(G, CS, state, fluxes) else ! This is needed because rigidity is potentially modified in the coupler. Reset ! in the ice shelf cavity: MJH + do j=isd,jed ; do i=isd,ied-1 ! changed stride fluxes%rigidity_ice_u(I,j) = (CS%kv_ice / CS%density_ice) * & min(CS%mass_shelf(i,j), CS%mass_shelf(i+1,j)) @@ -940,6 +945,10 @@ subroutine add_shelf_flux(G, CS, state, fluxes) endif if (associated(state%tauy_shelf)) then call vchksum(state%tauy_shelf, "tauy_shelf", G%HI, haloshift=0) + call vchksum(fluxes%rigidity_ice_u, "rigidity_ice_u", G%HI, haloshift=0) + call vchksum(fluxes%rigidity_ice_v, "rigidity_ice_v", G%HI, haloshift=0) + call vchksum(fluxes%frac_shelf_u, "frac_shelf_u", G%HI, haloshift=0) + call vchksum(fluxes%frac_shelf_v, "frac_shelf_v", G%HI, haloshift=0) endif endif @@ -967,7 +976,12 @@ subroutine add_shelf_flux(G, CS, state, fluxes) if ((asv1 + asv2 > 0.0) .and. associated(state%tauy_shelf)) & tauy2 = (asv1 * state%tauy_shelf(i,j-1)**2 + & asv2 * state%tauy_shelf(i,j)**2 ) / (asv1 + asv2) - fluxes%ustar(i,j) = MAX(CS%ustar_bg, sqrt(Irho0 * sqrt(taux2 + tauy2))) + ! GM, ustar must be the same as ustar_shelf, otherwise the results will + ! be different when starting from a RESTART file. + !fluxes%ustar(i,j) = MAX(CS%ustar_bg, sqrt(Irho0 * sqrt(taux2 + tauy2))) + + fluxes%ustar(i,j) = fluxes%ustar_shelf(i,j) + if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0 if (associated(fluxes%lw)) fluxes%lw(i,j) = 0.0 if (associated(fluxes%latent)) fluxes%latent(i,j) = 0.0 @@ -980,6 +994,7 @@ subroutine add_shelf_flux(G, CS, state, fluxes) endif endif + ! Add frazil formation diagnosed by the ocean model (J m-2) in the ! form of surface layer evaporation (kg m-2 s-1). Update lprec in the ! control structure for diagnostic purposes. @@ -1593,12 +1608,14 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti call register_restart_field(CS%taub_beta_eff_bilinear, vd, .true., CS%restart_CSp) endif - if (.not. solo_ice_sheet) then - vd = var_desc("ustar_shelf","m s-1","Friction velocity under ice shelves",z_grid='1') - call register_restart_field(fluxes%ustar_shelf, vd, .true., CS%restart_CSp) - vd = var_desc("iceshelf_melt","m year-1","Ice Shelf Melt Rate",z_grid='1') - call register_restart_field(fluxes%iceshelf_melt, vd, .true., CS%restart_CSp) - endif +! GM - I think we do not need to save ustar_shelf and iceshelf_melt in the restart file + +! if (.not. solo_ice_sheet) then +! vd = var_desc("ustar_shelf","m s-1","Friction velocity under ice shelves",z_grid='1') +! call register_restart_field(fluxes%ustar_shelf, vd, .true., CS%restart_CSp) +! vd = var_desc("iceshelf_melt","m year-1","Ice Shelf Melt Rate",z_grid='1') +! call register_restart_field(fluxes%iceshelf_melt, vd, .true., CS%restart_CSp) +! endif CS%restart_output_dir = dirs%restart_output_dir @@ -1648,8 +1665,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, & G, CS%restart_CSp) - - ! i think this call isnt necessary - all it does is set hmask to 3 at ! the dirichlet boundary, and now this is done elsewhere ! call initialize_shelf_mass(G, param_file, CS, .false.) From 0ea99d4595633d1886b879965d35c0d64db8cdd0 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 20 Sep 2016 16:22:12 -0400 Subject: [PATCH 41/46] Removed fluxes%ustar from add_shelf_flux fluxes%ustar is not used to compute melting. Instead, we use fluxes%ustar_shelf. Therefore, we do not need to update fluxes%ustar in the ice shelf code. --- src/ice_shelf/MOM_ice_shelf.F90 | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 6ffc8cde90..b5f6e6e03a 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -976,11 +976,11 @@ subroutine add_shelf_flux(G, CS, state, fluxes) if ((asv1 + asv2 > 0.0) .and. associated(state%tauy_shelf)) & tauy2 = (asv1 * state%tauy_shelf(i,j-1)**2 + & asv2 * state%tauy_shelf(i,j)**2 ) / (asv1 + asv2) - ! GM, ustar must be the same as ustar_shelf, otherwise the results will - ! be different when starting from a RESTART file. + + ! GM: melting is computed using ustar_shelf (and not ustar), which has already + ! been passed, so believe we do not need to update fluxes%ustar. !fluxes%ustar(i,j) = MAX(CS%ustar_bg, sqrt(Irho0 * sqrt(taux2 + tauy2))) - fluxes%ustar(i,j) = fluxes%ustar_shelf(i,j) if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0 if (associated(fluxes%lw)) fluxes%lw(i,j) = 0.0 @@ -1015,9 +1015,6 @@ subroutine add_shelf_flux(G, CS, state, fluxes) endif enddo ; enddo - if (CS%debug) then - call hchksum(fluxes%ustar_shelf, "ustar_shelf", G%HI, haloshift=0) - endif ! If the shelf mass is changing, the fluxes%rigidity_ice_[uv] needs to be ! updated here. From b004dc532a11c4c552d18e45b7f861a05d086541 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 21 Sep 2016 14:58:05 -0400 Subject: [PATCH 42/46] Many improvements in the ice shelf module * zero the fields of the ice shelf structure (ice_shelf_CS) that are used to compute melt. This is done right before melting is computed. This just affects diagnostics; * add safety checks; * pass fluxes%iceshelf_melt to the "melt" diagnostic. --- src/ice_shelf/MOM_ice_shelf.F90 | 47 ++++++++++++++++++++++----------- 1 file changed, 32 insertions(+), 15 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index b5f6e6e03a..7ef1b1befd 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -510,6 +510,14 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) iDens = 1.0/CS%density_ocean_avg + ! GM, zero some fields of the ice shelf structure (ice_shelf_CS) + ! these fields are already set to zero during initialization + ! However, they seem to be changed somewhere and, for diagnostic + ! reasons, it is better to set them to zero again. + CS%tflux_shelf(:,:) = 0.0; CS%exch_vel_t(:,:) = 0.0 + CS%lprec(:,:) = 0.0; CS%exch_vel_s(:,:) = 0.0 + CS%salt_flux(:,:) = 0.0; CS%t_flux(:,:) = 0.0 + CS%tfreeze(:,:) = 0.0 if (CS%shelf_mass_is_dynamic .and. CS%override_shelf_movement) & call update_shelf_mass(CS, Time) @@ -761,30 +769,38 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) enddo ! i-loop enddo ! j-loop + ! CS%lprec = precipitating liquid water into the ocean ( kg/(m^2 s) ) + ! melt in m/year + !fluxes%iceshelf_melt = CS%lprec * (86400.0*365.0/CS%density_ice) + ! Use rhoefw for ISOMIP experiments + fluxes%iceshelf_melt = CS%lprec * (86400.0*365.0/rho_fw) + do j=js,je do i=is,ie - ! GM: set lprec (therefore melting) to zero above a cutoff pressure - ! (CS%Rho0*CS%cutoff_depth*CS%g_Earth) this is needed for the isomip - ! test case. + ! GM: set melt to zero above a cutoff pressure + ! (CS%Rho0*CS%cutoff_depth*CS%g_Earth) this is needed for the isomip + ! test case. if ((CS%g_Earth * CS%mass_shelf(i,j)) < CS%Rho0*CS%cutoff_depth* & - CS%g_Earth) CS%lprec(i,j) = 0.0 - ! avoid negative fluxes - ! CS%lprec(i,j) = MAX(CS%lprec(i,j),0.0) - + CS%g_Earth) fluxes%iceshelf_melt(i,j) = 0.0 + + !GM, safety check + if (abs(fluxes%iceshelf_melt(i,j))>0.0) then + if (fluxes%ustar_shelf(i,j) == 0.0) then + write(*,*)'Something is wrong at i,j',i,j + call MOM_error(FATAL, & + "shelf_calc_flux: |melt| > 0 and star_shelf = 0.") + endif + endif enddo ! i-loop enddo ! j-loop - ! CS%lprec = precipitating liquid water into the ocean ( kg/(m^2 s) ) - ! melt in m/year - !fluxes%iceshelf_melt = CS%lprec * (86400.0*365.0/CS%density_ice) - ! Use rho_fw for ISOMIP experiments - fluxes%iceshelf_melt = CS%lprec * (86400.0*365.0/rho_fw) ! mass flux (kg/s), part of ISOMIP diags. ! first, allocate mass_flux. This is probably not the best way to do this. !im,jm=SHAPE(CS%lprec) - ALLOCATE ( mass_flux(G%ied,G%jed) ) - mass_flux = (CS%lprec) * CS%area_shelf_h + ALLOCATE ( mass_flux(G%ied,G%jed) ); mass_flux(:,:) = 0.0 + mass_flux = (CS%lprec) * CS%area_shelf_h + if (CS%DEBUG) then call hchksum (fluxes%iceshelf_melt, "melt rate", G%HI, haloshift=0) call hchksum (fluxes%ustar_shelf, "ustar_shelf calc", G%HI, haloshift=0) @@ -848,7 +864,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) if (CS%id_area_shelf_h > 0) call post_data(CS%id_area_shelf_h, CS%area_shelf_h, CS%diag) if (CS%id_ustar_shelf > 0) call post_data(CS%id_ustar_shelf, fluxes%ustar_shelf, CS%diag) ! Use freshwater density fotr ISOMIP exps. - if (CS%id_melt > 0) call post_data(CS%id_melt, (CS%lprec) * (86400.0*365.0/rho_fw), CS%diag) + !if (CS%id_melt > 0) call post_data(CS%id_melt, (CS%lprec) * (86400.0*365.0/rho_fw), CS%diag) + if (CS%id_melt > 0) call post_data(CS%id_melt, fluxes%iceshelf_melt, CS%diag) !if (CS%id_melt > 0) call post_data(CS%id_melt, (CS%lprec) * (86400.0*365.0/CS%density_ice), CS%diag) if (CS%id_thermal_driving > 0) call post_data(CS%id_thermal_driving, (state%sst-CS%tfreeze), CS%diag) if (CS%id_haline_driving > 0) call post_data(CS%id_haline_driving, (state%sss - Sbdry), CS%diag) From a12a67cc10a00515592d4f1934c00a03b39f2416 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 22 Sep 2016 10:43:29 -0400 Subject: [PATCH 43/46] Fixed ice shelf diag_to_Z - In the previous fix, the depth of the ice shelf (shelf_depth) was kept negative and it should be positive to dilate the water column properly. - No answer changes. --- src/diagnostics/MOM_diag_to_Z.F90 | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/diagnostics/MOM_diag_to_Z.F90 b/src/diagnostics/MOM_diag_to_Z.F90 index 74ea777154..40c3487134 100644 --- a/src/diagnostics/MOM_diag_to_Z.F90 +++ b/src/diagnostics/MOM_diag_to_Z.F90 @@ -232,7 +232,6 @@ subroutine calculate_Z_diag_fields(u, v, h, ssh_in, frac_shelf_h, dt, G, GV, CS) "diagnostic_fields_zstar: Module must be initialized before it is used.") ice_shelf = associated(frac_shelf_h) - shelf_depth(:) = 0. ! If no fields are needed, return if ((CS%id_u_z <= 0) .and. (CS%id_v_z <= 0) .and. (CS%num_tr_used < 1)) return @@ -246,15 +245,14 @@ subroutine calculate_Z_diag_fields(u, v, h, ssh_in, frac_shelf_h, dt, G, GV, CS) do j=js,je + shelf_depth(:) = 0. ! initially all is open ocean ! Remove all massless layers. do I=Isq,Ieq nk_valid(I) = 0 D_pt(I) = 0.5*(G%bathyT(i+1,j)+G%bathyT(i,j)) if (ice_shelf) then if (frac_shelf_h(i,j)+frac_shelf_h(i+1,j) > 0.) then ! under shelf - shelf_depth(I) = 0.5*(ssh(i+1,j)+ssh(i,j)) - else ! open ocean - shelf_depth(I) = 0.0 + shelf_depth(I) = abs(0.5*(ssh(i+1,j)+ssh(i,j))) endif endif enddo @@ -345,14 +343,13 @@ subroutine calculate_Z_diag_fields(u, v, h, ssh_in, frac_shelf_h, dt, G, GV, CS) enddo ; enddo ; enddo do J=Jsq,Jeq + shelf_depth(:) = 0.0 ! initially all is open ocean ! Remove all massless layers. do i=is,ie nk_valid(i) = 0 ; D_pt(i) = 0.5*(G%bathyT(i,j)+G%bathyT(i,j+1)) if (ice_shelf) then if (frac_shelf_h(i,j)+frac_shelf_h(i,j+1) > 0.) then ! under shelf - shelf_depth(i) = 0.5*(ssh(i,j)+ssh(i,j+1)) - else ! open ocean - shelf_depth(i) = 0.0 + shelf_depth(i) = abs(0.5*(ssh(i,j)+ssh(i,j+1))) endif endif enddo @@ -438,14 +435,13 @@ subroutine calculate_Z_diag_fields(u, v, h, ssh_in, frac_shelf_h, dt, G, GV, CS) enddo ; enddo ; enddo ; enddo do j=js,je + shelf_depth(:) = 0.0 ! initially all is open ocean ! Remove all massless layers. do i=is,ie nk_valid(i) = 0 ; D_pt(i) = G%bathyT(i,j) if (ice_shelf) then if (frac_shelf_h(i,j) > 0.) then ! under shelf - shelf_depth(i) = ssh(i,j) - else ! open ocean - shelf_depth(i) = 0.0 + shelf_depth(i) = abs(ssh(i,j)) endif endif enddo From c2115846d018eb4938b2a1afd3ded7fb76929682 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 13 Oct 2016 11:54:00 -0400 Subject: [PATCH 44/46] Removed the option to use usehmix when doing the surface averaged. --- src/core/MOM.F90 | 22 +--------------------- 1 file changed, 1 insertion(+), 21 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9419ebb3f5..f9d6b35da6 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -221,9 +221,6 @@ module MOM real :: Hmix_UV !< Depth scale over which to average surface flow to !! feedback to the coupler/driver (m). !! bulk mixed layer is not used. - logical :: usehmix !< If true, surface properties (e.g., SST) are averaged - !! over a distance given by HMIX_SFC_PROP. This parameter - !! is valid just when BULKMIXLAYER is defined. real :: missing=-1.0e34 !< missing data value for masked fields ! Flags needed to reach between start and finish phases of initialization @@ -1541,23 +1538,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) "case is the largest integer multiple of the coupling \n"//& "timestep that is less than or equal to DT_THERM.", default=.false.) - if (CS%bulkmixedlayer) then - call get_param(param_file, "MOM", "USE_HMIX_SFC_PROP", CS%usehmix, & - "Applied only if BULKMIXLAYER is used. If TRUE, surface properties \n"//& - "like SST, SSS, density and velocities are averaged over \n"//& - "a distance given by HMIX_SFC_PROP. These averaged fields are \n"//& - "passed to the surface state structure, which is then used \n"//& - "to, for example, compute melting under an ice shelf.", default=.false.) - if(CS%usehmix) then - call get_param(param_file, "MOM", "HMIX_SFC_PROP", CS%Hmix, & - "If USE_HMIX_SFC_PROP is true, HMIX_SFC_PROP is the depth \n"//& - "over which to average to find surface properties like \n"//& - "SST, SSS, density are surface velocities.", & - units="m", default=1.0) - endif - - endif - if (.not.CS%bulkmixedlayer) then call get_param(param_file, "MOM", "HMIX_SFC_PROP", CS%Hmix, & "If BULKMIXEDLAYER is false, HMIX_SFC_PROP is the depth \n"//& @@ -2933,7 +2913,7 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, GV, CS, p_atm) enddo ; enddo endif ; endif - if (CS%bulkmixedlayer .and. .not. CS%usehmix) then + if (CS%bulkmixedlayer) then state%SST => CS%tv%T(:,:,1) state%SSS => CS%tv%S(:,:,1) state%u => u(:,:,1) From 46ce1b03e0d83bb387d8b7ee8d84a0bc74122bc3 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 21 Oct 2016 11:30:49 -0400 Subject: [PATCH 45/46] Restored local work array needed for ISOMIP_tracers - For some unknown reason the declaration of a necessary work array was deleted. This code did not compile. --- src/tracer/ISOMIP_tracer.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index 9da0cdbb0c..a02be8368a 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -336,6 +336,7 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G real :: mmax real :: b1(SZI_(G)) ! b1 and c1 are variables used by the real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified real :: melt(SZI_(G),SZJ_(G)) ! melt water (positive for melting ! negative for freezing) integer :: i, j, k, is, ie, js, je, nz, m From 67f7c22bdd67c711c7892880d0935ae1750b3fcf Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 21 Oct 2016 11:36:35 -0400 Subject: [PATCH 46/46] Restored default value of S_REF to 35.0 - The default for S_REF in initialize_temp_salt_fit() was changed for no good reason (on branch gustavo-marques-dev/isomip_tmp) which broke several experiments that use this mode of initialization. --- src/initialization/MOM_state_initialization.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 16c2b198a6..8e4f014830 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1290,7 +1290,7 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, param_file, eqn_of_state, P_Ref units="degC", fail_if_missing=.true.) call get_param(param_file, mod, "S_REF", S_Ref, & "A reference salinity used in initialization.", units="PSU", & - default=33.8) + default=35.0) call get_param(param_file, mod, "FIT_SALINITY", fit_salin, & "If true, accept the prescribed temperature and fit the \n"//& "salinity; otherwise take salinity and fit temperature.", &