MODULE limsbc_2 !!====================================================================== !! *** MODULE limsbc_2 *** !! computation of the flux at the sea ice/ocean interface !!====================================================================== !! History : 1.0 ! 2000-01 (H. Goosse) Original code !! 2.0 ! 2002-07 (C. Ethe, G. Madec) re-writing F90 !! - ! 2006-07 (G. Madec) surface module !! - ! 2008-07 (C. Talandier,G. Madec) 2D fields for soce and sice !! 2.1 ! 2010-05 (Y. Aksenov G. Madec) salt flux + heat associated with emp !!---------------------------------------------------------------------- #if defined key_lim2 !!---------------------------------------------------------------------- !! 'key_lim2' LIM 2.0 sea-ice model !!---------------------------------------------------------------------- !! lim_sbc_2 : flux at the ice / ocean interface !!---------------------------------------------------------------------- USE par_oce ! ocean parameters USE dom_oce ! ocean domain USE sbc_ice ! ice surface boundary condition USE sbc_oce ! ocean surface boundary condition USE phycst ! physical constants USE albedo ! albedo parameters USE ice_2 ! LIM-2 sea-ice variables USE lbclnk ! ocean lateral boundary condition USE in_out_manager ! I/O manager USE iom ! USE prtctl ! Print control USE diaar5, ONLY : lk_diaar5 USE cpl_oasis3, ONLY : lk_cpl IMPLICIT NONE PRIVATE PUBLIC lim_sbc_2 ! called by sbc_ice_lim_2 REAL(wp) :: epsi16 = 1.e-16 ! constant values REAL(wp) :: rzero = 0.e0 REAL(wp) :: rone = 1.e0 REAL(wp), DIMENSION(jpi,jpj) :: soce_r, sice_r ! ocean and ice 2D constant salinity fields (used if lk_vvl=F) !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/LIM 3.3, UCL-LOCEAN-IPSL (2010) !! $Id$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_sbc_2( kt ) !!------------------------------------------------------------------- !! *** ROUTINE lim_sbc_2 *** !! !! ** Purpose : Update surface ocean boundary condition over areas !! that are at least partially covered by sea-ice !! !! ** Action : - comput. of the momentum, heat and freshwater/salt !! fluxes at the ice-ocean interface. !! - Update !! !! ** Outputs : - qsr : solar heat flux !! - qns : non-solar heat flux including heat content of mass flux !! - emp : mass flux !! - emps : salt flux due to Freezing/Melting !! - utau : sea surface i-stress (ocean referential) !! - vtau : sea surface j-stress (ocean referential) !! - fr_i : ice fraction !! - tn_ice : sea-ice surface temperature !! - alb_ice : sea-ice alberdo (lk_cpl=T) !! !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. !! Tartinville et al. 2001 Ocean Modelling, 3, 95-108. !!--------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! number of iteration !! INTEGER :: ji, jj ! dummy loop indices INTEGER :: ifvt, idfr , iadv, i1mfr ! local integers INTEGER :: iflt, ifrdv, ial , ifral ! - - INTEGER :: ii0, ii1, ij0, ij1 ! - - REAL(wp) :: zqsr, zqns, zqhc, zemp ! local scalars REAL(wp) :: zinda, zswitch, zcd ! - - REAL(wp) :: zfrldu, zutau, zu_io ! - - REAL(wp) :: zfrldv, zvtau, zv_io ! - - REAL(wp) :: zemp_snw, zfmm, zfsalt ! - - REAL(wp) :: zsang, zmod, zztmp, zfm REAL(wp), DIMENSION(jpi,jpj,1) :: zalb, zalbp ! 3D workspace REAL(wp), DIMENSION(jpi,jpj) :: ztio_u, ztio_v ! 2D workspace REAL(wp), DIMENSION(jpi,jpj) :: ztiomi, zqnsoce ! - - !!--------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition' IF(lwp) WRITE(numout,*) '~~~~~~~~~ ' ! ! 2D fields for constant ice and ocean salinities soce_r(:,:) = soce ! in order to use different value in the Baltic sea sice_r(:,:) = sice ! which is much less salty than polar regions ! IF( cp_cfg == "orca" ) THEN ! ORCA configuration IF( jp_cfg == 2 ) THEN ! ORCA_R2 configuration ii0 = 145 ; ii1 = 180 ! Baltic Sea ij0 = 113 ; ij1 = 130 ; soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 sice_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 !!gm add here the R1 R05 and R025 cases !! ELSEIF( jp_cfg == 1 ) THEN ! ORCA_R1 configuration !! ELSEIF( jp_cfg == 05 ) THEN ! ORCA_R05 configuration !! ELSEIF( jp_cfg == 025 ) THEN ! ORCA_R025 configuration !! !!gm or better introduce the baltic change as a function of lat/lon of the baltic sea !!end gm ENDIF ENDIF ! ENDIF zqnsoce(:,:) = qns(:,:) ! save non-solar flux prior to its modification by ice-ocean fluxes (diag.) ! zswitch = 1 ! standard levitating sea-ice : salt exchange only ! !!gm ice embedment ! SELECT CASE( nn_ice_embd ) ! levitating/embedded sea-ice ! CASE( 0 ) ; zswitch = 1 ! standard levitating sea-ice : salt exchange only ! CASE( 1, 2 ) ; zswitch = 0 ! other levitating sea-ice or embedded sea-ice : salt and volume fluxes ! END SELECT ! !!gm end embedment DO jj = 1, jpj DO ji = 1, jpi ! !------------------------------------------! ! ! heat flux at the ocean surface ! ! !------------------------------------------! zinda = 1.0 - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) ) ) ) ifvt = zinda * MAX( rzero , SIGN( rone, - phicif(ji,jj) ) ) i1mfr = 1.0 - MAX( rzero , SIGN( rone, - ( 1.0 - frld(ji,jj) ) ) ) idfr = 1.0 - MAX( rzero , SIGN( rone, frld(ji,jj) - pfrld(ji,jj) ) ) iflt = zinda * (1 - i1mfr) * (1 - ifvt ) ial = ifvt * i1mfr + ( 1 - ifvt ) * idfr iadv = ( 1 - i1mfr ) * zinda ifral = ( 1 - i1mfr * ( 1 - ial ) ) ifrdv = ( 1 - ifral * ( 1 - ial ) ) * iadv !!gm attempt to understand and comment the tricky flags used here.... ! !gm zinda = 1.0 - AINT( pfrld(ji,jj) ) ! = 0. free-ice ocean else 1. (after ice adv, but before ice thermo) !gm i1mfr = 1.0 - AINT( frld(ji,jj) ) ! = 0. free-ice ocean else 1. (after ice thermo) ! !gm IF( phicif(ji,jj) <= 0. ) THEN ; ifvt = zinda ! = 1. if (snow and no ice at previous time) else 0. ??? !gm ELSE ; ifvt = 0. ! correspond to a overmelting of snow in surface ablation !gm ENDIF ! ! !gm IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN ; idfr = 0. ! = 0. if lead fraction increases due to ice thermo !gm ELSE ; idfr = 1. !gm ENDIF ! !!$ iflt = zinda * (1 - i1mfr) * (1 - ifvt ) ! = 1. if ice (not only snow) at previous and pure ocean at current ! !!$ ial = ifvt * i1mfr + ( 1 - ifvt ) * idfr !!$! snow no ice ice ice or nothing lead fraction increases !!$! at previous now at previous !!$! -> ice aera increases ??? -> ice aera decreases ??? ! !!$ iadv = ( 1 - i1mfr ) * zinda !!$! pure ocean ice at !!$! at current previous !!$! -> = 1. if ice disapear between previous and current ! !!$ ifral = ( 1 - i1mfr * ( 1 - ial ) ) !!$! ice at ??? !!$! current !!$! -> ??? ! !!$ ifrdv = ( 1 - ifral * ( 1 - ial ) ) * iadv !!$! ice disapear ! ! ! - computation the solar flux at ocean surface #if defined key_coupled zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) ) #else zqsr = pfrld(ji,jj) * qsr(ji,jj) + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj) #endif ! ! - computation the non solar heat flux at ocean surface zqns = - ( 1. - thcm(ji,jj) ) * zqsr & ! part of the solar energy used in leads & + iflt * ( fscmbq(ji,jj) + ffltbif(ji,jj) ) & & + ifral * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdt_ice & & + ifrdv * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * r1_rdt_ice ! - store residual heat flux (put in the ocean at the next time-step) fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj) ! ??? ! ! - heat content of mass exchanged between ocean and sea-ice zqhc = ( rdq_snw(ji,jj) + rdq_ice(ji,jj) ) * r1_rdt_ice ! heat flux due to sown & ice heat content exchanges ! qsr(ji,jj) = zqsr ! solar heat flux qns(ji,jj) = zqns - fdtcn(ji,jj) + zqhc ! non solar heat flux ! !------------------------------------------! ! ! mass flux at the ocean surface ! ! !------------------------------------------! ! ! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area) #if defined key_coupled ! ! coupled mode: zemp = + emp_tot(ji,jj) & ! net mass flux over the grid cell (ice+ocean area) & - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) ) ! minus the mass flux intercepted by sea-ice #else ! ! forced mode: zemp = + emp(ji,jj) * frld(ji,jj) & ! mass flux over open ocean fraction & - tprecip(ji,jj) * ( 1. - frld(ji,jj) ) & ! liquid precip. over ice reaches directly the ocean & + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) ) & ! snow is intercepted by sea-ice (previous frld) #endif ! ! mass flux at the ocean/ice interface (sea ice fraction) zemp_snw = rdm_snw(ji,jj) * r1_rdt_ice ! snow melting = pure water that enters the ocean zfmm = rdm_ice(ji,jj) * r1_rdt_ice ! Freezing minus Melting (F-M) ! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s] zfsalt = - sice_r(ji,jj) * zfmm ! F-M salt exchange zcd = soce_r(ji,jj) * zfmm ! concentration/dilution term due to F-M ! ! salt flux only : add concentration dilution term in salt flux and no F-M term in volume flux ! salt and mass fluxes : non concentartion dilution term in salt flux and add F-M term in volume flux emps(ji,jj) = zfsalt + zswitch * zcd ! salt flux (+ C/D if no ice/ocean mass exchange) emp (ji,jj) = zemp + zemp_snw + ( 1.- zswitch) * zfmm ! mass flux (- F/M mass flux if no ice/ocean mass exchange) ! END DO END DO CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) ) CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) ) CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) ) IF( lk_diaar5 ) THEN CALL iom_put( 'isnwmlt_cea' , rdm_snw(:,:) * zrdtir ) CALL iom_put( 'fsal_virt_cea', soce_r(:,:) * rdm_ice(:,:) * zrdtir ) CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdm_ice(:,:) * zrdtir ) ENDIF !------------------------------------------! ! momentum flux at the ocean surface ! !------------------------------------------! IF ( ln_limdyn ) THEN ! Update the stress over ice-over area (only in ice-dynamic case) ! ! otherwise the atmosphere-ocean stress is used everywhere ! ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point) !CDIR NOVERRCHK DO jj = 1, jpj !CDIR NOVERRCHK DO ji = 1, jpi ! ... change the cosinus angle sign in the south hemisphere zsang = SIGN(1.e0, gphif(ji,jj) ) * sangvg ! ... ice velocity relative to the ocean at I-point zu_io = u_ice(ji,jj) - u_oce(ji,jj) zv_io = v_ice(ji,jj) - v_oce(ji,jj) zmod = SQRT( zu_io * zu_io + zv_io * zv_io ) zztmp = rhoco * zmod ! ... components of ice stress over ocean with a ice-ocean rotation angle (at I-point) ztio_u(ji,jj) = zztmp * ( cangvg * zu_io - zsang * zv_io ) ztio_v(ji,jj) = zztmp * ( cangvg * zv_io + zsang * zu_io ) ! ... module of ice stress over ocean (at I-point) ztiomi(ji,jj) = zztmp * zmod ! END DO END DO DO jj = 2, jpjm1 DO ji = 2, jpim1 ! NO vector opt. ! ... components of ice-ocean stress at U and V-points (from I-point values) zutau = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) ) zvtau = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) ) ! ... open-ocean (lead) fraction at U- & V-points (from T-point values) zfrldu = 0.5 * ( frld(ji,jj) + frld(ji+1,jj ) ) zfrldv = 0.5 * ( frld(ji,jj) + frld(ji ,jj+1) ) ! ... update components of surface ocean stress (ice-cover wheighted) utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * zutau vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * zvtau ! ... module of ice-ocean stress at T-points (from I-point values) zztmp = 0.25 * ( ztiomi(ji,jj) + ztiomi(ji+1,jj) + ztiomi(ji,jj+1) + ztiomi(ji+1,jj+1) ) ! ... update module of surface ocean stress (ice-cover wheighted) taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1. - frld(ji,jj) ) * zztmp ! END DO END DO CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary conditions CALL lbc_lnk( taum, 'T', 1. ) ! ENDIF !-----------------------------------------------! ! Storing the transmitted variables ! !-----------------------------------------------! !!gm where this is done ????? ==>>> limthd_2 not logic ??? !!gm fr_i(:,:) = 1.0 - frld(:,:) ! sea-ice fraction !!gm IF ( lk_cpl ) THEN ! coupled mode : tn_ice(:,:,1) = sist(:,:) ! sea-ice surface temperature ! ! snow/ice and ocean albedo CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb ) alb_ice(:,:,1) = 0.5 * ( zalbp(:,:,1) + zalb (:,:,1) ) ! Ice albedo (mean clear and overcast skys) ! CALL iom_put( "icealb_cea", alb_ice(:,:,1) * fr_i(:,:) ) ! ice albedo ENDIF IF(ln_ctl) THEN CALL prt_ctl(tab2d_1=qsr , clinfo1=' lim_sbc: qsr : ', tab2d_2=qns , clinfo2=' qns : ') CALL prt_ctl(tab2d_1=emp , clinfo1=' lim_sbc: emp : ', tab2d_2=emps , clinfo2=' emps : ') CALL prt_ctl(tab2d_1=utau , clinfo1=' lim_sbc: utau : ', mask1=umask, & & tab2d_2=vtau , clinfo2=' vtau : ' , mask2=vmask ) CALL prt_ctl(tab2d_1=fr_i , clinfo1=' lim_sbc: fr_i : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice : ') ENDIF ! END SUBROUTINE lim_sbc_2 #else !!---------------------------------------------------------------------- !! Default option : Dummy module NO LIM 2.0 sea-ice model !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_sbc_2 ! Dummy routine END SUBROUTINE lim_sbc_2 #endif !!====================================================================== END MODULE limsbc_2