MODULE iceupdate !!====================================================================== !! *** MODULE iceupdate *** !! Sea-ice : computation of the flux at the sea ice/ocean interface !!====================================================================== !! History : - ! 2006-07 (M. Vancoppelle) LIM3 original code !! 3.0 ! 2008-03 (C. Tallandier) surface module !! - ! 2008-04 (C. Tallandier) split in 2 + new ice-ocean coupling !! 3.3 ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea !! ! + simplification of the ice-ocean stress calculation !! 3.4 ! 2011-02 (G. Madec) dynamical allocation !! - ! 2012 (D. Iovino) salt flux change !! - ! 2012-05 (C. Rousset) add penetration solar flux !! 3.5 ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass !!---------------------------------------------------------------------- #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' ESIM sea-ice model !!---------------------------------------------------------------------- !! ice_update_alloc : allocate the iceupdate arrays !! ice_update_init : initialisation !! ice_update_flx : updates mass, heat and salt fluxes at the ocean surface !! ice_update_tau : update i- and j-stresses, and its modulus at the ocean surface !!---------------------------------------------------------------------- USE oce , ONLY : sshn, sshb USE phycst ! physical constants USE dom_oce ! ocean domain USE ice ! sea-ice: variables !!gm It should be probably better to pass these variable in argument of the routine, !!gm rather than having this long list in USE. This will also highlight what is updated, and what is just use. USE sbc_ice , ONLY : emp_oce, qns_oce, qsr_oce , qemp_oce , & & emp_ice, qsr_ice, qemp_ice, qevap_ice, alb_ice, tn_ice, cldf_ice, & & snwice_mass, snwice_mass_b, snwice_fmass USE sbc_oce , ONLY : nn_fsbc, ln_ice_embd, sfx, fr_i, qsr_tot, qns, qsr, fmmflx, emp, taum, utau, vtau !!gm end USE sbccpl ! Surface boundary condition: coupled interface USE icealb ! sea-ice: albedo parameters USE traqsr ! add penetration of solar flux in the calculation of heat budget USE icectl ! sea-ice: control prints USE bdy_oce , ONLY : ln_bdy ! USE in_out_manager ! I/O manager USE iom ! I/O manager library USE lib_mpp ! MPP library USE lib_fortran ! fortran utilities (glob_sum + no signed zero) USE lbclnk ! lateral boundary conditions (or mpp links) USE timing ! Timing IMPLICIT NONE PRIVATE PUBLIC ice_update_init ! called by ice_init PUBLIC ice_update_flx ! called by ice_stp PUBLIC ice_update_tau ! called by ice_stp REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_oce, vtau_oce ! air-ocean surface i- & j-stress [N/m2] REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmod_io ! modulus of the ice-ocean velocity [m/s] !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2017) !! $Id: iceupdate.F90 8411 2017-08-07 16:09:12Z clem $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION ice_update_alloc() !!------------------------------------------------------------------- !! *** ROUTINE ice_update_alloc *** !!------------------------------------------------------------------- ALLOCATE( utau_oce(jpi,jpj), vtau_oce(jpi,jpj), tmod_io(jpi,jpj), STAT=ice_update_alloc ) ! IF( lk_mpp ) CALL mpp_sum( ice_update_alloc ) IF( ice_update_alloc /= 0 ) CALL ctl_warn('ice_update_alloc: failed to allocate arrays') ! END FUNCTION ice_update_alloc SUBROUTINE ice_update_flx( kt ) !!------------------------------------------------------------------- !! *** ROUTINE ice_update_flx *** !! !! ** Purpose : Update the surface ocean boundary condition for heat !! salt and mass over areas where sea-ice is non-zero !! !! ** Action : - computes the heat and freshwater/salt fluxes !! at the ice-ocean interface. !! - Update the ocean sbc !! !! ** Outputs : - qsr : sea heat flux: solar !! - qns : sea heat flux: non solar !! - emp : freshwater budget: volume flux !! - sfx : salt flux !! - fr_i : ice fraction !! - tn_ice : sea-ice surface temperature !! - alb_ice : sea-ice albedo (recomputed only for coupled mode) !! !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. !! Tartinville et al. 2001 Ocean Modelling, 3, 95-108. !! These refs are now obsolete since everything has been revised !! The ref should be Rousset et al., 2015 !!--------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! number of iteration ! INTEGER :: ji, jj, jl, jk ! dummy loop indices REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2) REAL(wp) :: zqsr ! New solar flux received by the ocean REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalb_cs, zalb_os ! 3D workspace !!--------------------------------------------------------------------- IF( nn_timing == 1 ) CALL timing_start('ice_update') IF( kt == nit000 .AND. lwp ) THEN WRITE(numout,*) WRITE(numout,*)'ice_update_flx: update fluxes (mass, salt and heat) at the ice-ocean interface' WRITE(numout,*)'~~~~~~~~~~~~~~' ENDIF ! --- case we bypass ice thermodynamics --- ! IF( .NOT. ln_icethd ) THEN ! we suppose ice is impermeable => ocean is isolated from atmosphere hfx_in (:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:) hfx_out (:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) ftr_ice (:,:,:) = 0._wp emp_ice (:,:) = 0._wp qemp_ice (:,:) = 0._wp qevap_ice(:,:,:) = 0._wp ENDIF DO jj = 1, jpj DO ji = 1, jpi ! Solar heat flux reaching the ocean = zqsr (W.m-2) !--------------------------------------------------- zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) ) ) ! Total heat flux reaching the ocean = hfx_out (W.m-2) !--------------------------------------------------- zqmass = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr ! Add the residual from heat diffusion equation and sublimation (W.m-2) !---------------------------------------------------------------------- hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) + & & ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) ) ! New qsr and qns used to compute the oceanic heat flux at the next time step !---------------------------------------------------------------------------- qsr(ji,jj) = zqsr qns(ji,jj) = hfx_out(ji,jj) - zqsr ! Mass flux at the atm. surface !----------------------------------- wfx_sub(ji,jj) = wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj) ! Mass flux at the ocean surface !------------------------------------ ! case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED) ! ------------------------------------------------------------------------------------- ! The idea of this approach is that the system that we consider is the ICE-OCEAN system ! Thus FW flux = External ( E-P+snow melt) ! Salt flux = Exchanges in the ice-ocean system then converted into FW ! Associated to Ice formation AND Ice melting ! Even if i see Ice melting as a FW and SALT flux ! ! mass flux from ice/ocean wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) & & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) IF ( ln_pnd_fw ) wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj) ! add the snow melt water to snow mass flux to the ocean wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj) ! mass flux at the ocean/ice interface fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) ) ! F/M mass flux save at least for biogeochemical model emp(ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! mass flux + F/M mass flux (always ice/ocean mass exchange) ! Salt flux at the ocean surface !------------------------------------------ sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj) & & + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj) ! Mass of snow and ice per unit area !---------------------------------------- snwice_mass_b(ji,jj) = snwice_mass(ji,jj) ! save mass from the previous ice time step ! ! new mass per unit area snwice_mass (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) ) ! ! time evolution of snow+ice mass snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice END DO END DO ! Storing the transmitted variables !---------------------------------- fr_i (:,:) = at_i(:,:) ! Sea-ice fraction tn_ice(:,:,:) = t_su(:,:,:) ! Ice surface temperature ! Snow/ice albedo (only if sent to coupler, useless in forced mode) !------------------------------------------------------------------ CALL ice_alb( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos ! alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) ! IF( lrst_ice ) THEN !* write snwice_mass fields in the restart file CALL update_rst( 'WRITE', kt ) ENDIF ! ! controls IF( ln_icediachk .AND. .NOT. ln_bdy) CALL ice_cons_final('iceupdate') ! conservation IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints IF( ln_ctl ) CALL ice_prt3D ('iceupdate') ! prints IF( nn_timing == 1 ) CALL timing_stop ('ice_update') ! timing ! END SUBROUTINE ice_update_flx SUBROUTINE ice_update_tau( kt, pu_oce, pv_oce ) !!------------------------------------------------------------------- !! *** ROUTINE ice_update_tau *** !! !! ** Purpose : Update the ocean surface stresses due to the ice !! !! ** Action : * at each ice time step (every nn_fsbc time step): !! - compute the modulus of ice-ocean relative velocity !! (*rho*Cd) at T-point (C-grid) or I-point (B-grid) !! tmod_io = rhoco * | U_ice-U_oce | !! - update the modulus of stress at ocean surface !! taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce | !! * at each ocean time step (every kt): !! compute linearized ice-ocean stresses as !! Utau = tmod_io * | U_ice - pU_oce | !! using instantaneous current ocean velocity (usually before) !! !! NB: - ice-ocean rotation angle no more allowed !! - here we make an approximation: taum is only computed every ice time step !! This avoids mutiple average to pass from T -> U,V grids and next from U,V grids !! to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton... !! !! ** Outputs : - utau, vtau : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes !! - taum : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes !!--------------------------------------------------------------------- INTEGER , INTENT(in) :: kt ! ocean time-step index REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zat_u, zutau_ice, zu_t, zmodt ! local scalar REAL(wp) :: zat_v, zvtau_ice, zv_t, zrhoco ! - - !!--------------------------------------------------------------------- IF( nn_timing == 1 ) CALL timing_start('ice_update') IF( kt == nit000 .AND. lwp ) THEN WRITE(numout,*) WRITE(numout,*)'ice_update_tau: update stress at the ice-ocean interface' WRITE(numout,*)'~~~~~~~~~~~~~~' ENDIF zrhoco = rau0 * rn_cio ! IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step) DO jj = 2, jpjm1 !* update the modulus of stress at ocean surface (T-point) DO ji = fs_2, fs_jpim1 ! ! 2*(U_ice-U_oce) at T-point zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj) zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) ! ! |U_ice-U_oce|^2 zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t ) ! ! update the ocean stress modulus taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * zrhoco * zmodt tmod_io(ji,jj) = zrhoco * SQRT( zmodt ) ! rhoco * |U_ice-U_oce| at T-point END DO END DO CALL lbc_lnk_multi( taum, 'T', 1., tmod_io, 'T', 1. ) ! utau_oce(:,:) = utau(:,:) !* save the air-ocean stresses at ice time-step vtau_oce(:,:) = vtau(:,:) ! ENDIF ! ! !== every ocean time-step ==! ! DO jj = 2, jpjm1 !* update the stress WITHOUT a ice-ocean rotation angle DO ji = fs_2, fs_jpim1 ! Vect. Opt. zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp ! ice area at u and V-points zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp ! ! linearized quadratic drag formulation zutau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) ) zvtau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) ) ! ! stresses at the ocean surface utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice END DO END DO CALL lbc_lnk_multi( utau, 'U', -1., vtau, 'V', -1. ) ! lateral boundary condition ! IF( nn_timing == 1 ) CALL timing_stop('ice_update') ! END SUBROUTINE ice_update_tau SUBROUTINE ice_update_init !!------------------------------------------------------------------- !! *** ROUTINE ice_update_init *** !! !! ** Purpose : ??? !! !!------------------------------------------------------------------- INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar !!------------------------------------------------------------------- ! IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'ice_update_init: ???? ' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' ! ! allocate ice_update array IF( ice_update_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_update_init : unable to allocate standard arrays' ) ! CALL update_rst( 'READ' ) !* read or initialize all required files ! END SUBROUTINE ice_update_init SUBROUTINE update_rst( cdrw, kt ) !!--------------------------------------------------------------------- !! *** ROUTINE rhg_evp_rst *** !! !! ** Purpose : Read or write RHG file in restart file !! !! ** Method : use of IOM library !!---------------------------------------------------------------------- CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step ! INTEGER :: iter ! local integer INTEGER :: id1 ! local integer !!---------------------------------------------------------------------- ! IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialize ! ! --------------- IF( ln_rstart ) THEN !* Read the restart file ! id1 = iom_varid( numrir, 'snwice_mass' , ldstop = .FALSE. ) ! IF( id1 > 0 ) THEN ! fields exist CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass' , snwice_mass ) CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b ) ELSE ! start from rest IF(lwp) WRITE(numout,*) ' ==>> previous run without snow-ice mass output then set it' snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) snwice_mass_b(:,:) = snwice_mass(:,:) ENDIF ELSE !* Start from rest IF(lwp) WRITE(numout,*) ' ==>> start from rest: set the snow-ice mass' snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) snwice_mass_b(:,:) = snwice_mass(:,:) ENDIF ! ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file ! ! ------------------- IF(lwp) WRITE(numout,*) '---- update-rst ----' iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1 ! CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass' , snwice_mass ) CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b ) ! ENDIF ! END SUBROUTINE update_rst #else !!---------------------------------------------------------------------- !! Default option Dummy module NO ESIM sea-ice model !!---------------------------------------------------------------------- #endif !!====================================================================== END MODULE iceupdate