hosing experiments subroutines

src/fesom_module.F90

@@ -534,7 +534,8 @@ end subroutine fesom_init
   subroutine fesom_runloop(current_nsteps)
     use fesom_main_storage_module
 !   use openacc_lib
-    integer, intent(in) :: current_nsteps 
+    integer, intent(in) :: current_nsteps
+    real(kind=WP)                            :: hSv
     ! EO parameters
     integer n, nstart, ntotal, tr_num
 
@@ -690,6 +691,12 @@ subroutine fesom_runloop(current_nsteps)
             if (flag_debug .and. f%mype==0)  print *, achar(27)//'[34m'//' --> call oce_fluxes_mom...'//achar(27)//'[0m'
             call oce_fluxes_mom(f%ice, f%dynamics, f%partit, f%mesh) ! momentum only
             call oce_fluxes(f%ice, f%dynamics, f%tracers, f%partit, f%mesh)
+
+            !___freshwater depth hosing routine_______________________________________
+            !
+            hSv=0.1 !define freshwater anomaly magnitude
+            call fw_depth_anomaly(f%tracers%data(2)%values, f%tracers%data(1)%values,  hSv, f%partit, f%mesh)
+
         end if
         call before_oce_step(f%dynamics, f%tracers, f%partit, f%mesh) ! prepare the things if required
         f%t2 = MPI_Wtime()

src/ice_oce_coupling.F90

@@ -43,6 +43,23 @@ subroutine oce_fluxes_mom(ice, dynamics, partit, mesh)
         type(t_partit), intent(inout), target :: partit
         type(t_mesh)  , intent(in)   , target :: mesh
         end subroutine oce_fluxes_mom
+
+        subroutine fw_surf_anomaly(hSv, tracers, partit, mesh)
+        use o_PARAM
+        use o_ARRAYS
+        USE MOD_TRACER
+        USE MOD_PARTIT
+        use MOD_MESH
+        use MOD_PARSUP
+        USE g_CONFIG
+        use g_comm_auto
+        use g_support
+        type(t_tracer), intent(in),    target :: tracers
+        type(t_partit), intent(inout), target :: partit
+        type(t_mesh)  , intent(in)   , target :: mesh
+        real(kind=WP) , intent(in)            :: hSv
+        end subroutine
+
     end interface
 end module oce_fluxes_interface
 
@@ -279,7 +296,7 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh)
     type(t_mesh)  , intent(in)   , target :: mesh
     !___________________________________________________________________________
     integer                    :: n, elem, elnodes(3),n1
-    real(kind=WP)              :: rsss, net
+    real(kind=WP)              :: rsss, net, hSv
     real(kind=WP), allocatable :: flux(:)
     !___________________________________________________________________________
     real(kind=WP), dimension(:,:), pointer :: temp, salt
@@ -398,6 +415,11 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh)
 !$OMP END PARALLEL DO
 #endif
 
+    !___freshwater surface hosing routine_______________________________________
+    !
+    hSv=0.1 !define freshwater anomaly magnitude
+    call fw_surf_anomaly(hSv, tracers, partit, mesh)
+
     if (use_icebergs) then
         call icb2fesom(mesh, partit, ice)
     end if
@@ -839,3 +861,215 @@ end subroutine oce_fluxes
 !
 !
 !_______________________________________________________________________________
+subroutine fw_surf_anomaly(hSv, tracers, partit, mesh)  !!!! subroutine to apply freshwater at the surface
+
+    use o_PARAM
+    use o_ARRAYS
+    USE MOD_TRACER
+    USE MOD_PARTIT
+    use MOD_MESH
+    use MOD_PARSUP
+    USE g_CONFIG
+    use g_comm_auto
+    use g_support
+    implicit none
+
+    type(t_partit), intent(inout), target  :: partit
+    type(t_mesh),   intent(in),    target  :: mesh
+    type(t_tracer), intent(in),    target  :: tracers
+    real(kind=WP),  intent(in)             :: hSv
+    real(kind=WP)                          :: x, y, net
+    integer                                :: n, ed(2)
+    real(kind=WP), dimension(:,:), pointer :: temp
+    real(kind=WP), parameter               :: lhf = 334000. !Latent heat of fusion for sea ice (J/Kg)
+    real(kind=WP), parameter               :: rhofwt = 1000. !Fresh water density (Kg/m3)
+#include "associate_part_def.h"
+#include "associate_mesh_def.h"
+#include "associate_part_ass.h"
+#include "associate_mesh_ass.h"
+      temp          => tracers%data(1)%values(:,:)
+
+      hosing_flux=0._WP
+      hosing_heat_flux=0._WP
+
+     !!freshwater anomaly is applied on all surface nodes south of 60S (uncomment the next 4 lines)
+     !do n=1, myDim_nod2D+eDim_nod2D
+     !   y = geo_coord_nod2D(2,n)/rad
+     !   if (y<-60._WP) hosing_flux(n)=1._WP
+     !end do
+
+     !!freshwater anomaly is applied only on surface coastal nodes south of 60S (uncomment the next 6 lines)
+     !!the mask is created looping on all the mesh edges and considering only the nodes belonging to non internal edges 
+     do n=1, myDim_edge2D
+        ed=mesh%edges(:, n)
+        if (myList_edge2D(n) <= mesh%edge2D_in) cycle
+        y = sum(geo_coord_nod2D(2,ed))/2._WP/rad
+        if (y<-60._WP) hosing_flux(ed)=1._WP
+     end do
+
+     !!this is to increase the number of adjacent nodes on which the anomaly is applied, starting from the coastal nodes. 
+     !!We use first a smoothing function, choosing the number of nodes (here 8) to smooth the freshwater anomaly and then equalize the value for all nodes
+     !call smooth_nod (hosing_flux, 8, partit, mesh) 
+     !do n=1,myDim_nod2d+eDim_nod2D
+     !   if (hosing_flux(n)>0.0_WP) hosing_flux(n)=1.0_WP
+     !end do
+
+     !!computing the area-integrated value of hosing_flux (net is the spatial integral) and then normalizing it so that its area integral equals 1
+     !!and then multiplying for the actual freshwater flux to impose (hSv)
+     call integrate_nod(hosing_flux, net, partit, mesh)
+
+     if (abs(net)>1.e-6) then
+        hosing_flux=hosing_flux/net*hSv*1.e6 !the freshwater anomaly is transformed from Sv in m/s (Sv*1.e6 =  m/s)
+     end if
+
+     water_flux=water_flux-hosing_flux !freshwater flux is negative going into the ocean, so the sign is negative 
+
+     !!Computing the heat used to melt the ice corresponding to the freshwater flux
+     !!mass_wat = hosing_flux*rhofwt !mass flux of fresh water (Kg/m2*s) 
+     !!the heat flux is in W/m2 so no multiplication for the node area is needed 
+     !!mass_ice = hosing_flux_ice*rhoice =  hosing_flux*rhofwt !mass of melted ice (Kg/m2*s)
+     !! in therms of mass they are equal
+     !!hosing_latent_heat_flux = mass_ice*Lf !(J/m2*s=W/m2) 
+
+     !!to remove latent heat of fusion in the freshwater experiment uncomment the next line
+     !hosing_heat_flux = hosing_flux*rhofwt*lhf !(W/m2)
+
+     !!Computing the heat flux change due to the freshwater temperature 
+     !!by default freshwater temperature is assumed equal to surrounding temperature so heat flux should be removed if we want the freshwater to have a different temperature
+     !!since we are adding super cool freshwater (water at freezing temp) we have a heat flux going out of the ocean, i.e. positive)   
+     !!with t_freezing=-1.8, (temp(1,:)-t_freezing) is always positive (temp is almost always >-1.8, so hosing_heat_flux>0)
+
+     !!to modify latent heat to take into account different freshwater temperature uncomment the next 2 lines
+     !hosing_heat_flux=(temp(1,:)+1.8)*hosing_flux*vcpw  
+     !heat_flux=heat_flux+hosing_heat_flux !(W/m2) !heat flux is negative going into the ocean, the hosing heat flux goes out of the ocean, so the sign is positive 
+
+!!print the freshwater amount in the log file (it is printed for every processor so it is a lot or printouts)     
+!write(*,*) mype, 'hSv=', hSv, 'net=', net, 'hf=', minval(hosing_flux, 1), maxval(hosing_flux, 1)
+
+end subroutine fw_surf_anomaly
+!
+!
+!_______________________________________________________________________________
+subroutine fw_depth_anomaly(tts, ttt, hSv, partit, mesh) !!!! subroutine to apply freshwater at a chosen depth
+    use o_PARAM
+    use o_ARRAYS
+    USE MOD_PARTIT
+    use MOD_MESH
+    use MOD_PARSUP
+    use g_comm_auto
+    use g_support
+    use o_tracers
+    use g_config,         only: dt
+    implicit none
+
+    integer                               :: n, ed(2), row, k, nz, nlev1, nlev2, nzmin, nzmax, elem1, elem2
+    real(kind=WP)                         :: x, y, net, spar(100)
+    real(kind=WP),  intent(in)            :: hSv
+    type(t_mesh),   intent(in),    target :: mesh
+    type(t_partit), intent(inout), target :: partit
+    real(kind=WP),  intent (inout)        :: tts(mesh%nl-1,partit%myDim_nod2D+partit%eDim_nod2D)
+    real(kind=WP),  intent (inout)        :: ttt(mesh%nl-1,partit%myDim_nod2D+partit%eDim_nod2D)
+    real(kind=WP), save, allocatable      :: mask(:), mask3D(:,:)
+    real(kind=WP), parameter              :: lhf = 334000. !Latent heat of fusion for sea ice (J/Kg)
+    real(kind=WP), parameter              :: rhofwt = 1000. !Fresh water density (Kg/m3)
+    logical, save                         :: lfirst=.true.
+#include "associate_part_def.h"
+#include "associate_mesh_def.h"
+#include "associate_part_ass.h"
+#include "associate_mesh_ass.h"
+    hosing_flux3D=0._WP
+    hosing_heat_flux3D=0._WP
+
+    !!generating the 3D mask to apply the anomaly (only for the first iteration)
+    if (lfirst) then
+       allocate(mask(myDim_nod2D+eDim_nod2D), mask3D(mesh%nl-1, myDim_nod2D+eDim_nod2D))
+       mask =0._WP
+       mask3D=0._WP
+       spar=0._WP
+       !!the 2D mask is created in the same way as for the surface application first selecting nodes belonging to coastal edges
+       !!and then using the smoothing function to increase it including 5 rows of adjacent nodes
+       do n=1, myDim_edge2D
+          ed=mesh%edges(:, n)
+          if (myList_edge2D(n) <= mesh%edge2D_in) cycle
+             y = sum(geo_coord_nod2D(2,ed))/2._WP/rad
+          if (y<-60._WP) mask(ed)=1._WP
+       end do
+
+       call smooth_nod (mask, 5, partit, mesh)
+       where (mask>1.e-12)
+             mask=1.0_WP
+       end where
+
+       !!create a 3D mask to apply the freshwater anomaly between the chosen levels nzmin and nzmax 
+       do n=1, myDim_edge2D
+          ed = mesh%edges(:, n)
+          !!select all the elements adjacent to the edges with both nodes in the 2D mask
+          if (any(mask(ed)<1.e-12)) cycle 
+          elem1 = mesh%edge_tri(1, n)
+          elem2 = mesh%edge_tri(2, n)
+          !!determine how many vertical levels exist for each element
+          nlev1 = mesh%nlevels(elem1)
+          nlev2 = 1
+          if (elem2 > 0) nlev2 = mesh%nlevels(elem2)
+          !!choose a vertical range of levels where the freshwater anomaly will be applied
+          nzmin = 16
+          nzmax = 19
+          !!select all nodes belonging to edges fully inside the mask2D in the coosen vertical levels
+          do k = nzmin, nzmax
+             if (k < nlev1 .or. k < nlev2) then !!Apply mask only to levels that exist in at least one of the two elements
+             mask3D(k, ed) = 1.0_WP
+             end if
+          end do
+
+       end do
+
+       !!collapse the 3D mask back to a 2D mask (a node is kept in mask only if it has at least one active depth level in mask3D)
+       !!to compute the area integral and normalize the mask horizontally
+       do row=1,myDim_nod2d+eDim_nod2D
+          if (sum(mask3D(:, row))<1.e-12) mask(row)=0._WP
+       end do
+
+       call integrate_nod(mask, net, partit, mesh)
+       if (abs(net)>1.e-6) then
+          mask=mask/net*1.e6 !hSv will be applied later as it changes in time
+       end if
+
+       !!Building the 3D distribution
+       do row=1,myDim_nod2d+eDim_nod2D 
+          spar(:mesh%nl)=0._WP
+          !!build a volume-weighted vertical shape inside the chosen levels (to distribute the anomaly proportionally to the volume)
+          do k=ulevels_nod2D(row), nlevels_nod2D(row)-1
+             spar(k)=areasvol(k,row)*hnode(k,row)*mask3D(k, row)                           
+          end do
+          !!normalize the vertical distribution so it sums to 1
+          spar(:nlevels_nod2D(row)-1)=spar(:nlevels_nod2D(row)-1)/max(sum(spar(:nlevels_nod2D(row)-1)), 1.e-12)
+          !!convert the vertically normalized weights so that ts vertical integral at each node matches the horizontal forcing defined by mask(row) 
+          do k=ulevels_nod2D(row), nlevels_nod2D(row)-1
+             spar(k)=spar(k)/areasvol(k,row)/hnode(k,row)*mask(row)*area(1,row)*dt
+          end do
+          mask3D(:nlevels_nod2D(row)-1,row)=spar(:nlevels_nod2D(row)-1)
+       end do
+
+       lfirst=.false. !set to false so that the mask generation is not repeated
+    !!printout the area/volume-integrated flux to check it is correct
+    !call integrate_nod(mask3D, net, partit, mesh)
+    !write(*,*) 'integrated mask3D=', net
+    !call integrate_nod(mask, net, partit, mesh)
+    !write(*,*) 'integrated mask2D=', net
+    end if
+
+    !!updating the salinity field 
+    !!(Snew = Sold (1-dS) where dS=mask3D*hSv is the fraction of the gridcell volume where saltwater is replaced by freshwater during the timestep)
+    tts=tts-mask3D*hSv*tts
+    !call integrate_nod(tts, net, partit, mesh)
+    !write(*,*) 'integrated salinity=', net
+
+    !hosing_heat_flux3D = mask3D*hsv*rhofwt*lhf    
+    !ttt = ttt-hosing_heat_flux3D/vcpw !temperature field modification to remove latent heat of fusion
+    !ttt=ttt-mask3D*hSv*(ttt+1.8) !temperature field modification to have meltwaterwater added at -1.8 deg.
+    !call integrate_nod(ttt, net, partit, mesh)
+    !write(*,*) 'integrated temperature=', net
+
+    hosing_flux3D=mask3D*hSv
+
+end subroutine fw_depth_anomaly

src/io_meandata.F90

@@ -617,6 +617,17 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh)
 CASE ('ice_rejectsalt')
     if (SPP) call def_stream(nod2D , myDim_nod2D , 'ice_rejectsalt' , 'salt flux from plum parameterisation ', 'm/s*psu', ice_rejected_salt(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
 
+!___________________________________________________________________________________________________________________________________
+! output hosing experiment  
+CASE ('hfw       ')
+call def_stream(nod2D, myDim_nod2D,     'hfw',      'hosing water flux',               'm/s',    hosing_flux(:),            io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
+CASE ('hfh       ')
+call def_stream(nod2D, myDim_nod2D,     'hfh',      'hosing heat flux',                'W',      hosing_heat_flux(:),       io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
+CASE ('hfw3D     ')
+    call def_stream((/nl-1, nod2D/),  (/nl-1, myDim_nod2D/),  'hfw3D',     'hosing water flux in 3D', 'm/s',    hosing_flux3D(:,:),          io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
+CASE ('hfh3D     ')
+    call def_stream((/nl-1, nod2D/),  (/nl-1, myDim_nod2D/),  'hfh3D',     'hosing heat flux in 3D',  'W',      hosing_heat_flux3D(:,:),     io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
+
 !___________________________________________________________________________________________________________________________________
 ! output KPP vertical mixing schemes
 CASE ('kpp_obldepth   ')

src/io_restart.F90

@@ -210,13 +210,13 @@ subroutine ini_ocean_io(dynamics, tracers, partit, mesh)
          write(longname,'(A15,i4.4)') 'passive tracer ', j
          units='none'
      END SELECT
-     if ((tracers%data(j)%ID==101) .or. (tracers%data(j)%ID==102) .or. (tracers%data(j)%ID==103)) then
+     if ((tracers%data(j)%ID==101) .or. (tracers%data(j)%ID==102) .or. (tracers%data(j)%ID==103) .or. (tracers%data(j)%ID==304)) then
         call oce_files%def_node_var_optional(trim(trname), trim(longname), trim(units), tracers%data(j)%values(:,:), mesh, partit)
      else
         call oce_files%def_node_var(trim(trname), trim(longname), trim(units), tracers%data(j)%values(:,:), mesh, partit)
      endif
      longname=trim(longname)//', Adams-Bashforth'
-     if ((tracers%data(j)%ID==101) .or. (tracers%data(j)%ID==102) .or. (tracers%data(j)%ID==103)) then
+     if ((tracers%data(j)%ID==101) .or. (tracers%data(j)%ID==102) .or. (tracers%data(j)%ID==103) .or. (tracers%data(j)%ID==304)) then
         call oce_files%def_node_var_optional(trim(trname)//'_AB', trim(longname), trim(units), tracers%data(j)%valuesAB(:,:),    mesh, partit)
      else
         call oce_files%def_node_var(trim(trname)//'_AB', trim(longname), trim(units), tracers%data(j)%valuesAB(:,:),    mesh, partit)

src/oce_ale_tracer.F90

@@ -1027,7 +1027,17 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh)
         !  (BUT CHECK!)              |    |                         |    |
         !                            v   (+)                        v   (+)
         !
-        tr(nzmin)= tr(nzmin)+bc_surface(n, tracers%data(tr_num)%ID, trarr(nzmin,n), nzmin, partit, mesh, sst(nzmin,n), sss(nzmin,n), a_ice(n))
+
+        !for depth freshwater exp
+        if (tracers%data(tr_num)%ID==304) then
+            do nz=nzmin,nzmax
+                tr(nz)= tr(nz)+bc_surface(n, tracers%data(tr_num)%ID, trarr(nz,n), nz, partit, mesh, sst(nz,n), sss(nz,n), a_ice(n))
+            end do
+        else
+            tr(nzmin)= tr(nzmin)+bc_surface(n, tracers%data(tr_num)%ID, trarr(nzmin,n), nzmin, partit, mesh, sst(nzmin,n), sss(nzmin,n), a_ice(n))
+        end if
+         !for surface freshwater exp 
+        tr(nzmin)= tr(nzmin)+bc_surface(n, tracers%data(tr_num)%ID, trarr(nzmin,n), nzmin, partit, mesh, sst(nzmin,n), sss(nzmin,n), a_ice(n)) 
 
         !_______________________________________________________________________
         ! The forward sweep algorithm to solve the three-diagonal matrix
@@ -1660,6 +1670,9 @@ FUNCTION bc_surface(n, id, sval, nzmin, partit, mesh, sst, sss, aice)
         bc_surface=0.0_WP
     CASE (303)
         bc_surface=0.0_WP
+    CASE (304)
+        bc_surface= dt*(hosing_flux3D(nzmin,n)) !for depth freshwater exp
+        bc_surface= dt*(hosing_flux(n)) !for surf freshwater exp
     CASE (501) ! ice-shelf water due to basal melting
         if (nzmin==1) then
            bc_surface = 0.0_WP

src/oce_modules.F90

@@ -225,6 +225,8 @@ MODULE o_ARRAYS
 real(kind=WP), allocatable         :: heat_flux_in(:) !to keep the unmodified (by SW penetration etc.) heat flux 
 real(kind=WP), allocatable         :: Tsurf_ib(:) ! kh 15.03.21 additional array for asynchronous iceberg computations
 real(kind=WP), allocatable    :: water_flux(:), fw_ice(:), fw_snw(:), Ssurf(:)
+real(kind=WP), allocatable    :: hosing_flux(:), hosing_heat_flux(:)
+real(kind=WP), allocatable    :: hosing_flux3D(:,:), hosing_heat_flux3D(:,:)
 real(kind=WP), allocatable    :: Ssurf_ib(:) ! kh 15.03.21 additional array for asynchronous iceberg computations
 real(kind=WP), allocatable    :: virtual_salt(:), relax_salt(:)
 real(kind=WP), allocatable    :: Tclim(:,:), Sclim(:,:)

src/oce_setup_step.F90

@@ -857,6 +857,8 @@ SUBROUTINE arrays_init(num_tracers, partit, mesh)
     allocate(relax2clim(node_size)) 
     allocate(heat_flux(node_size), Tsurf(node_size))
     allocate(water_flux(node_size), Ssurf(node_size))
+    allocate(hosing_flux(node_size), hosing_heat_flux(node_size))
+    allocate(hosing_flux3D(nl-1,node_size), hosing_heat_flux3D(nl-1,node_size))
     allocate(fw_ice(node_size), fw_snw(node_size))
     allocate(relax_salt(node_size))
     allocate(virtual_salt(node_size))
@@ -952,6 +954,10 @@ SUBROUTINE arrays_init(num_tracers, partit, mesh)
     Tsurf=0.0_WP
 
     water_flux=0.0_WP
+    hosing_flux=0.0_WP
+    hosing_heat_flux=0.0_WP
+    hosing_flux3D=0.0_WP
+    hosing_heat_flux3D=0.0_WP
     fw_ice    =0.0_WP
     fw_snw    =0.0_WP
     relax_salt=0.0_WP
@@ -1377,7 +1383,16 @@ SUBROUTINE oce_initial_state(tracers, partit, mesh)
                 write (id_string, "(I3)") id
                 write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string)
             end if
-
+
+        !_______________________________________________________________________
+        CASE (304) ! passive tracer for water hosing experiment
+            tracers%data(i)%values(:,:)=0.0_WP
+            if (mype==0) then
+                write (i_string,  "(I3)") i
+                write (id_string, "(I3)") id
+                write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string)
+            end if
+
         !_______________________________________________________________________
         CASE (501) ! ice-shelf water due to basal melting
             tracers%data(i)%values(:,:)=0.0_WP