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' LIM 3.0 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 par_oce ! ocean parameters 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 ! albedo parameters USE traqsr ! add penetration of solar flux in the calculation of heat budget USE domvvl ! Variable volume USE icectl ! ??? USE bdy_oce , ONLY : ln_bdy ! USE in_out_manager ! I/O manager USE iom ! xIO server USE lbclnk ! ocean lateral boundary condition - MPP exchanges USE lib_mpp ! MPP library USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 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] REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: soce_0 , sice_0 ! cst SSS and ice salinity (levitating sea-ice) !! * 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( soce_0(jpi,jpj) , utau_oce(jpi,jpj) , & & sice_0(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_flx') IF( kt == nit000 .AND. lwp ) THEN WRITE(numout,*) WRITE(numout,*)'ice_update_flx' WRITE(numout,*)'~~~~~~~~~~~~~~' ENDIF ! --- case we bypass ice thermodynamics --- ! IF( .NOT. ln_limthd ) 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(:,:,:) ! ! conservation test IF( ln_limdiachk .AND. .NOT. ln_bdy) CALL ice_cons_final( 'iceupdate' ) ! ! control prints IF( ln_limctl ) CALL ice_prt( kt, iiceprt, jiceprt, 3, ' - Final state ice_update - ' ) IF( ln_ctl ) CALL ice_prt3D( 'iceupdate' ) ! IF( nn_timing == 1 ) CALL timing_stop('ice_update_flx') ! 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_tau') IF( kt == nit000 .AND. lwp ) THEN WRITE(numout,*) WRITE(numout,*)'ice_update_tau' 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_tau') ! 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 : sea-ice ????' 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' ) ! soce_0(:,:) = soce ! constant SSS and ice salinity used in levitating case 0 (i.e. virtual salt flux) sice_0(:,:) = sice WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND. & ! reduced values in the Baltic Sea area & 54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp ) soce_0(:,:) = 4._wp sice_0(:,:) = 2._wp END WHERE ! IF( .NOT.ln_rstart ) THEN ! set ! snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) ! snow+ice mass snwice_mass_b(:,:) = snwice_mass(:,:) ! IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 !!gm I really don't like this stuff here... Find a way to put that elsewhere or differently !!gm IF( .NOT.ln_linssh ) THEN DO jk = 1,jpkm1 ! adjust initial vertical scale factors e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) ) ) e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshb(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) ) ) END DO e3t_a(:,:,:) = e3t_b(:,:,:) !!gm we are in no-restart case, so sshn=sshb ==>> faster calculation: !! REAL(wp) :: ze3t ! local scalar !! REAL(wp), DIMENSION(jpi,jpj) :: z2d ! workspace !! !! WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:) !! ELSEWHERE ; z2d(:,:) = 1._wp !! END WHERE !! DO jk = 1,jpkm1 ! adjust initial vertical scale factors !! ze3t = e3t_0(:,:,jk) * z2d(:,:) !! e3t_n(:,:,jk) = ze3t !! e3t_b(:,:,jk) = ze3t !! e3t_a(:,:,jk) = ze3t !! END DO !!gm but since it is only done at the initialisation.... just the following can be acceptable: ! DO jk = 1,jpkm1 ! adjust initial vertical scale factors ! e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1)) ) ! END DO ! e3t_b(:,:,:) = e3t_n(:,:,:) ! e3t_a(:,:,:) = e3t_n(:,:,:) !!gm end ! Reconstruction of all vertical scale factors at now and before time-steps ! ========================================================================= ! Horizontal scale factor interpolations ! -------------------------------------- CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) ! Vertical scale factor interpolations ! ------------------------------------ CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) ! t- and w- points depth ! ---------------------- !!gm not sure of that.... gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) gdepw_n(:,:,1) = 0.0_wp gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) DO jk = 2, jpk gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk ) gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1) gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:) END DO ENDIF ENDIF ENDIF ! .NOT. ln_rstart ! END SUBROUTINE ice_update_init #else !!---------------------------------------------------------------------- !! Default option Dummy module NO LIM3 sea-ice model !!---------------------------------------------------------------------- #endif !!====================================================================== END MODULE iceupdate