Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • nemo/nemo
  • sparonuz/nemo
  • hatfield/nemo
  • extdevs/nemo
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 962 additions and 742 deletions
src/ABL/sbcabl.F90
View file @ 3efbcb17
......@@ -44,6 +44,8 @@ MODULE sbcabl
PUBLIC sbc_abl_init ! routine called in sbcmod module
PUBLIC sbc_abl ! routine called in sbcmod module
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OPA 3.7 , NEMO-consortium (2014)
!! $Id: sbcabl.F90 6416 2016-04-01 12:22:17Z clem $
......@@ -308,8 +310,8 @@ CONTAINS
!! - Perform 1 time-step of the ABL model
!! - Finalize flux computation in blk_oce_2
!!
!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
!! - vtau : j-component of the stress at V-point (N/m2)
!! ** Outputs : - utau : i-component of the stress at T-point (N/m2)
!! - vtau : j-component of the stress at T-point (N/m2)
!! - taum : Wind stress module at T-point (N/m2)
!! - wndm : Wind speed module at T-point (m/s)
!! - qsr : Solar heat flux over the ocean (W/m2)
......@@ -340,7 +342,7 @@ CONTAINS
CALL blk_oce_1( kt, u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in
& tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), & ! <<= in
& sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m , & ! <<= in
& sf(jp_slp )%fnow(:,:,1) , sst_m(A2D(0)), ssu_m(A2D(0)), ssv_m(A2D(0)), & ! <<= in
& sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in
& sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) , & ! <<= in
& tsk_m, zssq, zcd_du, zsen, zlat, zevp ) ! =>> out
......@@ -348,7 +350,7 @@ CONTAINS
#if defined key_si3
CALL blk_ice_1( u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in
& tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), & ! <<= in
& sf(jp_slp)%fnow(:,:,1) , u_ice, v_ice, tm_su , & ! <<= in
& sf(jp_slp)%fnow(:,:,1) , u_ice(A2D(0)), v_ice(A2D(0)), tm_su , & ! <<= in
& pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui ) ! <<= out
#endif
......
src/ICE/ice.F90
View file @ 3efbcb17
......@@ -451,6 +451,8 @@ MODULE ice
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qcn_ice_bot !: Bottom conduction flux (W/m2)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qcn_ice_top !: Surface conduction flux (W/m2)
!
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: ice.F90 15388 2021-10-17 11:33:47Z clem $
......@@ -464,71 +466,103 @@ CONTAINS
!!-----------------------------------------------------------------
INTEGER :: ice_alloc
!
INTEGER :: ierr(16), ii
INTEGER :: ierr(21), ii
!!-----------------------------------------------------------------
ierr(:) = 0
ii = 0
! ----------------- !
! == FULL ARRAYS == !
! ----------------- !
! * Ice global state variables
ii = ii + 1
ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , STAT=ierr(ii) )
ii = 1
ALLOCATE( u_oce (jpi,jpj) , v_oce (jpi,jpj) , ht_i_new (jpi,jpj) , fraz_frac (jpi,jpj) , &
& strength (jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) , &
& delta_i (jpi,jpj) , divu_i (jpi,jpj) , shear_i (jpi,jpj) , &
& aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv (jpi,jpj) , STAT=ierr(ii) )
ii = ii + 1
ALLOCATE( h_i (jpi,jpj,jpl) , a_i (jpi,jpj,jpl) , v_i (jpi,jpj,jpl) , &
& v_s (jpi,jpj,jpl) , h_s (jpi,jpj,jpl) , &
& s_i (jpi,jpj,jpl) , sv_i(jpi,jpj,jpl) , o_i (jpi,jpj,jpl) , oa_i (jpi,jpj,jpl) , &
& a_ip (jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl), &
& v_il (jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , &
& t_su (jpi,jpj,jpl) , t_s (jpi,jpj,nlay_s,jpl) , t_i(jpi,jpj,nlay_i,jpl) , sz_i(jpi,jpj,nlay_i,jpl) , &
& ato_i(jpi,jpj) , STAT = ierr(ii) )
ii = ii + 1
ALLOCATE( t_bo (jpi,jpj) , wfx_snw_sni(jpi,jpj) , &
& wfx_snw (jpi,jpj) , wfx_snw_dyn(jpi,jpj) , wfx_snw_sum(jpi,jpj) , wfx_snw_sub(jpi,jpj) , &
& wfx_ice (jpi,jpj) , wfx_sub (jpi,jpj) , wfx_ice_sub(jpi,jpj) , wfx_lam (jpi,jpj) , &
& wfx_pnd (jpi,jpj) , &
& wfx_bog (jpi,jpj) , wfx_dyn (jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) , &
& wfx_res (jpi,jpj) , wfx_sni (jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) , &
& rn_amax_2d (jpi,jpj) , &
& qsb_ice_bot(jpi,jpj) , qlead (jpi,jpj) , &
& sfx_res (jpi,jpj) , sfx_bri (jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) , &
& sfx_bog (jpi,jpj) , sfx_bom (jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) , &
& hfx_res (jpi,jpj) , hfx_snw (jpi,jpj) , hfx_sub(jpi,jpj) , &
& qt_atm_oi (jpi,jpj) , qt_oce_ai (jpi,jpj) , fhld (jpi,jpj) , &
& hfx_sum (jpi,jpj) , hfx_bom (jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , &
& hfx_opw (jpi,jpj) , hfx_thd (jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) , &
& hfx_err_dif(jpi,jpj) , wfx_err_sub(jpi,jpj) , STAT=ierr(ii) )
ALLOCATE( e_s(jpi,jpj,nlay_s,jpl) , e_i(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) )
! * Ice global state variables
! * Before values of global variables
ii = ii + 1
ALLOCATE( qtr_ice_bot(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) , &
& h_i (jpi,jpj,jpl) , a_i (jpi,jpj,jpl) , v_i (jpi,jpj,jpl) , &
& v_s (jpi,jpj,jpl) , h_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) , &
& s_i (jpi,jpj,jpl) , sv_i (jpi,jpj,jpl) , o_i (jpi,jpj,jpl) , &
& oa_i (jpi,jpj,jpl) , bv_i (jpi,jpj,jpl) , STAT=ierr(ii) )
ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) )
! * fluxes
ii = ii + 1
ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , &
& vt_i (jpi,jpj) , vt_s (jpi,jpj) , st_i(jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) , &
& et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s(jpi,jpj) , &
& sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s(jpi,jpj) , &
& om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj), icb_mask(jpi,jpj), STAT=ierr(ii) )
ALLOCATE( wfx_res(jpi,jpj) , sfx_res(jpi,jpj) , hfx_res(jpi,jpj) , STAT=ierr(ii) ) ! full arrays since it is used in conjonction with global variables
! * ice rheology
ii = ii+1
ALLOCATE( u_oce (jpi,jpj) , v_oce (jpi,jpj) , &
& strength (jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) , &
& aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv (jpi,jpj) , &
& icb_mask (jpi,jpj) , STAT=ierr(ii) )
! * mean and total
ii = ii + 1
ALLOCATE( vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , & ! full arrays since they are used in rheology
& vt_ip(jpi,jpj) , vt_il(jpi,jpj) , at_ip(jpi,jpj) , STAT=ierr(ii) )
! * others
ii = ii + 1
ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) )
ALLOCATE( t_bo(jpi,jpj) , rn_amax_2d(jpi,jpj) , STAT=ierr(ii) )
! -------------------- !
! == REDUCED ARRAYS == !
! -------------------- !
! * Ice global state variables
ii = ii + 1
ALLOCATE( t_i(jpi,jpj,nlay_i,jpl) , e_i(jpi,jpj,nlay_i,jpl) , sz_i(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) )
ALLOCATE( bv_i(A2D(0),jpl) , a_ip_frac(A2D(0),jpl) , a_ip_eff(A2D(0),jpl) , STAT=ierr(ii) )
! * Before values of global variables
ii = ii + 1
ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl), &
& v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , &
& dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) , STAT = ierr(ii) )
ALLOCATE( at_i_b(A2D(0)) , h_i_b (A2D(0),jpl) , a_i_b(A2D(0),jpl) , v_i_b(A2D(0),jpl) , &
& v_s_b (A2D(0),jpl) , h_s_b (A2D(0),jpl) , &
& v_ip_b(A2D(0),jpl) , v_il_b(A2D(0),jpl) , &
& sv_i_b(A2D(0),jpl) , e_i_b (A2D(0),nlay_i,jpl) , e_s_b(A2D(0),nlay_s,jpl) , STAT=ierr(ii) )
! * fluxes
ii = ii + 1
ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , hm_il(jpi,jpj) , vt_il(jpi,jpj) , STAT = ierr(ii) )
ALLOCATE( qsb_ice_bot(A2D(0)) , qlead (A2D(0)) , qt_atm_oi (A2D(0)) , qt_oce_ai (A2D(0)) , fhld (A2D(0)) , &
& wfx_snw_sni(A2D(0)) , wfx_snw (A2D(0)) , wfx_snw_dyn(A2D(0)) , wfx_snw_sum(A2D(0)) , wfx_snw_sub(A2D(0)) , &
& wfx_ice (A2D(0)) , wfx_sub (A2D(0)) , wfx_ice_sub(A2D(0)) , wfx_lam (A2D(0)) , &
& wfx_pnd (A2D(0)) , &
& wfx_bog (A2D(0)) , wfx_dyn (A2D(0)) , wfx_bom(A2D(0)) , wfx_sum(A2D(0)) , &
& wfx_sni (A2D(0)) , wfx_opw (A2D(0)) , wfx_spr(A2D(0)) , &
& &
& sfx_bri (A2D(0)) , sfx_dyn (A2D(0)) , sfx_sub(A2D(0)) , sfx_lam(A2D(0)) , &
& sfx_bog (A2D(0)) , sfx_bom (A2D(0)) , sfx_sum(A2D(0)) , sfx_sni(A2D(0)) , sfx_opw(A2D(0)) , &
& hfx_snw (A2D(0)) , hfx_sub (A2D(0)) , &
& hfx_sum (A2D(0)) , hfx_bom (A2D(0)) , hfx_bog(A2D(0)) , hfx_dif(A2D(0)) , &
& hfx_opw (A2D(0)) , hfx_thd (A2D(0)) , hfx_dyn(A2D(0)) , hfx_spr(A2D(0)) , &
& hfx_err_dif(A2D(0)) , wfx_err_sub(A2D(0)) , STAT=ierr(ii) )
ii = ii + 1
ALLOCATE( qtr_ice_bot(A2D(0),jpl) , cnd_ice(A2D(0),jpl) , t1_ice(A2D(0),jpl) , STAT=ierr(ii) )
! * ice rheology
ii = ii+1
ALLOCATE( delta_i(A2D(0)) , divu_i(A2D(0)) , shear_i(A2D(0)) , STAT=ierr(ii) )
! * Old values of global variables
! * mean and total
ii = ii + 1
ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl) , h_i_b(jpi,jpj,jpl), &
& v_ip_b(jpi,jpj,jpl) , v_il_b(jpi,jpj,jpl) , &
& a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , &
& STAT=ierr(ii) )
ALLOCATE( st_i (A2D(0)) , et_i (A2D(0)) , et_s(A2D(0)) , hm_i (A2D(0)) , &
& hm_ip(A2D(0)) , hm_il(A2D(0)) , tm_i(A2D(0)) , tm_s (A2D(0)) , &
& sm_i (A2D(0)) , hm_s (A2D(0)) , om_i(A2D(0)) , bvm_i(A2D(0)) , &
& tm_su(A2D(0)) , STAT=ierr(ii) )
! * others
ii = ii + 1
ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , at_i_b(jpi,jpj) , STAT=ierr(ii) )
ALLOCATE( tau_icebfr(A2D(0)) , dh_i_sum_2d(A2D(0),jpl) , dh_s_mlt_2d(A2D(0),jpl) , STAT=ierr(ii) )
ii = 1
ALLOCATE( ht_i_new (A2D(0)) , fraz_frac (A2D(0)) , STAT=ierr(ii) )
! * Ice thickness distribution variables
ii = ii + 1
......@@ -536,19 +570,19 @@ CONTAINS
! * Ice diagnostics
ii = ii + 1
ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), &
& diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), &
& diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), diag_vpnd(jpi,jpj), &
& diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) )
ALLOCATE( diag_trp_vi (A2D(0)) , diag_trp_vs (A2D(0)) , diag_trp_ei (A2D(0)) , &
& diag_trp_es (A2D(0)) , diag_trp_sv (A2D(0)) , diag_heat (A2D(0)) , &
& diag_sice (A2D(0)) , diag_vice (A2D(0)) , diag_vsnw (A2D(0)) , diag_aice(A2D(0)) , diag_vpnd(A2D(0)), &
& diag_adv_mass(A2D(0)) , diag_adv_salt(A2D(0)) , diag_adv_heat(A2D(0)) , STAT=ierr(ii) )
! * Ice conservation
ii = ii + 1
ALLOCATE( diag_v (jpi,jpj) , diag_s (jpi,jpj) , diag_t (jpi,jpj), &
& diag_fv(jpi,jpj) , diag_fs(jpi,jpj) , diag_ft(jpi,jpj), STAT=ierr(ii) )
ALLOCATE( diag_v (A2D(0)) , diag_s (A2D(0)) , diag_t (A2D(0)), &
& diag_fv(A2D(0)) , diag_fs(A2D(0)) , diag_ft(A2D(0)), STAT=ierr(ii) )
! * SIMIP diagnostics
ii = ii + 1
ALLOCATE( t_si(jpi,jpj,jpl) , tm_si(jpi,jpj) , qcn_ice_bot(jpi,jpj,jpl) , qcn_ice_top(jpi,jpj,jpl) , STAT = ierr(ii) )
ALLOCATE( t_si(A2D(0),jpl) , tm_si(A2D(0)) , qcn_ice_bot(A2D(0),jpl) , qcn_ice_top(A2D(0),jpl) , STAT = ierr(ii) )
ice_alloc = MAXVAL( ierr(:) )
IF( ice_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_alloc: failed to allocate arrays.' )
......
src/ICE/icealb.F90
View file @ 3efbcb17
......@@ -48,7 +48,7 @@ MODULE icealb
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, pcloud_fra, palb_ice )
SUBROUTINE ice_alb( ld_pnd_alb, pt_su, ph_ice, ph_snw, pafrac_pnd, ph_pnd, pcloud_fra, palb_ice )
!!----------------------------------------------------------------------
!! *** ROUTINE ice_alb ***
!!
......@@ -94,16 +94,16 @@ CONTAINS
!! Brandt et al. 2005, J. Climate, vol 18
!! Grenfell & Perovich 2004, JGR, vol 109
!!----------------------------------------------------------------------
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: pt_su ! ice surface temperature (Kelvin)
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: ph_ice ! sea-ice thickness
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: ph_snw ! snow depth
LOGICAL , INTENT(in ) :: ld_pnd_alb ! effect of melt ponds on albedo
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: pafrac_pnd ! melt pond relative fraction (per unit ice area)
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: ph_pnd ! melt pond depth
REAL(wp), INTENT(in ), DIMENSION(:,:) :: pcloud_fra ! cloud fraction
REAL(wp), INTENT( out), DIMENSION(:,:,:) :: palb_ice ! albedo of ice
LOGICAL , INTENT(in ) :: ld_pnd_alb ! effect of melt ponds on albedo
REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: pt_su ! ice surface temperature (Kelvin)
REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: ph_ice ! sea-ice thickness
REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: ph_snw ! snow depth
REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: pafrac_pnd ! melt pond relative fraction (per unit ice area)
REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: ph_pnd ! melt pond depth
REAL(wp), INTENT(in ), DIMENSION(A2D(0)) :: pcloud_fra ! cloud fraction
REAL(wp), INTENT( out), DIMENSION(A2D(0),jpl) :: palb_ice ! albedo of ice
!
REAL(wp), DIMENSION(jpi,jpj,jpl) :: za_s_fra ! ice fraction covered by snow
REAL(wp), DIMENSION(A2D(0),jpl) :: za_s_fra ! ice fraction covered by snow
INTEGER :: ji, jj, jl ! dummy loop indices
REAL(wp) :: z1_c1, z1_c2,z1_c3, z1_c4 ! local scalar
REAL(wp) :: z1_href_pnd ! inverse of the characteristic length scale (Lecomte et al. 2015)
......@@ -121,10 +121,10 @@ CONTAINS
z1_c3 = 1._wp / 0.02_wp
z1_c4 = 1._wp / 0.03_wp
!
CALL ice_var_snwfra( ph_snw, za_s_fra ) ! calculate ice fraction covered by snow
CALL ice_var_snwfra( ph_snw(:,:,:), za_s_fra(:,:,:) ) ! calculate ice fraction covered by snow
!
DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! palb_ice used over the full domain in icesbc
DO_2D( 0, 0, 0, 0 ) ! palb_ice used over the full domain in icesbc
!
!---------------------------------------------!
!--- Specific snow, ice and pond fractions ---!
......@@ -164,10 +164,10 @@ CONTAINS
zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd(ji,jj,jl) * z1_href_pnd )
!
! !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions
zalb_os = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1)
zalb_os = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * smask0(ji,jj)
!
zalb_cs = zalb_os - ( - 0.1010_wp * zalb_os * zalb_os &
& + 0.1933_wp * zalb_os - 0.0148_wp ) * tmask(ji,jj,1)
& + 0.1933_wp * zalb_os - 0.0148_wp ) * smask0(ji,jj)
!
! albedo depends on cloud fraction because of non-linear spectral effects
palb_ice(ji,jj,jl) = ( 1._wp - pcloud_fra(ji,jj) ) * zalb_cs + pcloud_fra(ji,jj) * zalb_os
......
src/ICE/icectl.F90
View file @ 3efbcb17
......@@ -83,25 +83,33 @@ CONTAINS
CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
REAL(wp) , INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
!!
INTEGER :: ji, jj, jl ! dummy loop index
REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat
REAL(wp), DIMENSION(jpi,jpj,10) :: ztmp3
REAL(wp), DIMENSION(jpi,jpj,jpl,8) :: ztmp4
REAL(wp), DIMENSION(10) :: zchk3
REAL(wp), DIMENSION(8) :: zchk4
REAL(wp), DIMENSION(A2D(0),10) :: ztmp3
REAL(wp), DIMENSION(A2D(0),jpl,8) :: ztmp4
REAL(wp), DIMENSION(10) :: zchk3
REAL(wp), DIMENSION(8) :: zchk4
!!-------------------------------------------------------------------
!
! -- quantities -- !
ztmp3(:,:,1) = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ! volume
ztmp3(:,:,2) = SUM( sv_i * rhoi, dim=3 ) * e1e2t ! salt
ztmp3(:,:,3) = ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ! heat
!
! -- fluxes -- !
ztmp3(:,:,4) = ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd & ! mass
& + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) * e1e2t
ztmp3(:,:,5) = ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw & ! salt
& + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t
ztmp3(:,:,6) = ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & ! heat
& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t
DO_2D( 0, 0, 0, 0 )
! -- quantities -- !
ztmp3(ji,jj,1) = SUM( v_i(ji,jj,:) * rhoi + v_s(ji,jj,:) * rhos + &
& ( v_ip(ji,jj,:) + v_il(ji,jj,:) ) * rhow ) * e1e2t(ji,jj) ! volume
ztmp3(ji,jj,2) = SUM( sv_i(ji,jj,:) * rhoi ) * e1e2t(ji,jj) ! salt
ztmp3(ji,jj,3) = ( SUM( SUM( e_i(ji,jj,:,:), dim=2 ) ) + & ! heat
& SUM( SUM( e_s(ji,jj,:,:), dim=2 ) ) ) * e1e2t(ji,jj)
!
! -- fluxes -- !
ztmp3(ji,jj,4) = ( wfx_bog (ji,jj) + wfx_bom (ji,jj) + wfx_sum (ji,jj) + wfx_sni (ji,jj) & ! mass
& + wfx_opw (ji,jj) + wfx_res (ji,jj) + wfx_dyn (ji,jj) + wfx_lam (ji,jj) + wfx_pnd(ji,jj) &
& + wfx_snw_sni(ji,jj) + wfx_snw_sum(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sub(ji,jj) &
& + wfx_ice_sub(ji,jj) + wfx_spr(ji,jj) ) * e1e2t(ji,jj)
ztmp3(ji,jj,5) = ( sfx_bri(ji,jj) + sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj) & ! salt
& + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj) ) * e1e2t(ji,jj)
ztmp3(ji,jj,6) = ( hfx_sum(ji,jj) + hfx_bom(ji,jj) + hfx_bog(ji,jj) + hfx_dif(ji,jj) + hfx_opw(ji,jj) + hfx_snw(ji,jj) & ! heat
& - hfx_thd(ji,jj) - hfx_dyn(ji,jj) - hfx_res(ji,jj) - hfx_sub(ji,jj) - hfx_spr(ji,jj) ) * e1e2t(ji,jj)
!
END_2D
!
! -- global sum -- !
zchk3(1:6) = glob_sum_vec( 'icectl', ztmp3(:,:,1:6) )
......@@ -123,25 +131,33 @@ CONTAINS
zdiag_heat = ( zchk3(3) - pdiag_t ) * r1_Dt_ice + ( zchk3(6) - pdiag_ft )
! -- max concentration diag -- !
ztmp3(:,:,7) = SUM( a_i, dim=3 )
zchk3(7) = glob_max( 'icectl', ztmp3(:,:,7) )
DO_2D( 0, 0, 0, 0 )
ztmp3(ji,jj,7) = SUM( a_i(ji,jj,:) )
END_2D
zchk3(7) = glob_max( 'icectl', ztmp3(:,:,7) )
! -- advection scheme is conservative? -- !
ztmp3(:,:,8 ) = diag_adv_mass * e1e2t
ztmp3(:,:,9 ) = diag_adv_heat * e1e2t
ztmp3(:,:,10) = SUM( a_i + epsi10, dim=3 ) * e1e2t ! ice area (+epsi10 to set a threshold > 0 when there is no ice)
zchk3(8:10) = glob_sum_vec( 'icectl', ztmp3(:,:,8:10) )
DO_2D( 0, 0, 0, 0 )
ztmp3(ji,jj,8 ) = diag_adv_mass(ji,jj) * e1e2t(ji,jj)
ztmp3(ji,jj,9 ) = diag_adv_heat(ji,jj) * e1e2t(ji,jj)
ztmp3(ji,jj,10) = SUM( a_i(ji,jj,:) + epsi10 ) * e1e2t(ji,jj) ! ice area (+epsi10 to set a threshold > 0 when there is no ice)
END_2D
zchk3(8:10) = glob_sum_vec( 'icectl', ztmp3(:,:,8:10) )
! -- min diags -- !
ztmp4(:,:,:,1) = v_i
ztmp4(:,:,:,2) = v_s
ztmp4(:,:,:,3) = v_ip
ztmp4(:,:,:,4) = v_il
ztmp4(:,:,:,5) = a_i
ztmp4(:,:,:,6) = sv_i
ztmp4(:,:,:,7) = SUM( e_i, dim=3 )
ztmp4(:,:,:,8) = SUM( e_s, dim=3 )
zchk4(1:8) = glob_min_vec( 'icectl', ztmp4(:,:,:,1:8) )
DO jl = 1, jpl
DO_2D( 0, 0, 0, 0 )
ztmp4(ji,jj,jl,1) = v_i(ji,jj,jl)
ztmp4(ji,jj,jl,2) = v_s(ji,jj,jl)
ztmp4(ji,jj,jl,3) = v_ip(ji,jj,jl)
ztmp4(ji,jj,jl,4) = v_il(ji,jj,jl)
ztmp4(ji,jj,jl,5) = a_i(ji,jj,jl)
ztmp4(ji,jj,jl,6) = sv_i(ji,jj,jl)
ztmp4(ji,jj,jl,7) = SUM( e_i(ji,jj,:,jl) )
ztmp4(ji,jj,jl,8) = SUM( e_s(ji,jj,:,jl) )
END_2D
ENDDO
zchk4(1:8) = glob_min_vec( 'icectl', ztmp4(:,:,:,1:8) )
IF( lwp ) THEN
! check conservation issues
......@@ -188,17 +204,21 @@ CONTAINS
!!-------------------------------------------------------------------
CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
!!
REAL(wp), DIMENSION(jpi,jpj,4) :: ztmp
REAL(wp), DIMENSION(4) :: zchk
INTEGER :: ji, jj ! dummy loop index
REAL(wp), DIMENSION(A2D(0),4) :: ztmp
REAL(wp), DIMENSION(4) :: zchk
!!-------------------------------------------------------------------
ztmp(:,:,1) = ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + wfx_pnd + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ! mass diag
ztmp(:,:,2) = ( sfx + diag_sice - diag_adv_salt ) * e1e2t ! salt
ztmp(:,:,3) = ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ! heat
! equivalent to this:
!! ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
!! & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t )
ztmp(:,:,4) = SUM( a_i + epsi10, dim=3 ) * e1e2t ! ice area (+epsi10 to set a threshold > 0 when there is no ice)
DO_2D( 0, 0, 0, 0 )
!
ztmp(ji,jj,1) = ( wfx_ice (ji,jj) + wfx_snw (ji,jj) + wfx_pnd (ji,jj) + wfx_spr(ji,jj) + wfx_sub(ji,jj) &
& + diag_vice(ji,jj) + diag_vsnw(ji,jj) + diag_vpnd(ji,jj) - diag_adv_mass(ji,jj) ) * e1e2t(ji,jj) ! mass diag
ztmp(ji,jj,2) = ( sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) ) * e1e2t(ji,jj) ! salt
ztmp(ji,jj,3) = ( qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj) ) * e1e2t(ji,jj) ! heat
! equivalent to this:
!! ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
!! & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t )
ztmp(ji,jj,4) = SUM( a_i(ji,jj,:) + epsi10 ) * e1e2t(ji,jj) ! ice area (+epsi10 to set a threshold > 0 when there is no ice)
END_2D
! global sums
zchk(1:4) = glob_sum_vec( 'icectl', ztmp(:,:,1:4) )
......@@ -226,11 +246,11 @@ CONTAINS
!!-------------------------------------------------------------------
INTEGER , INTENT(in) :: icount ! called at: =0 the begining of the routine, =1 the end
CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
REAL(wp) , DIMENSION(jpi,jpj), INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
REAL(wp) , DIMENSION(A2D(0)), INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
!!
REAL(wp), DIMENSION(jpi,jpj) :: zdiag_mass, zdiag_salt, zdiag_heat, &
& zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax
INTEGER :: jl, jk
REAL(wp), DIMENSION(A2D(0)) :: zdiag_mass, zdiag_salt, zdiag_heat, &
& zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax
INTEGER :: ji, jj, jl, jk
LOGICAL :: ll_stop_m = .FALSE.
LOGICAL :: ll_stop_s = .FALSE.
LOGICAL :: ll_stop_t = .FALSE.
......@@ -239,62 +259,79 @@ CONTAINS
!
IF( icount == 0 ) THEN
pdiag_v = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 )
pdiag_s = SUM( sv_i * rhoi , dim=3 )
pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 )
DO_2D( 0, 0, 0, 0 )
pdiag_v(ji,jj) = SUM( v_i(ji,jj,:) * rhoi + v_s(ji,jj,:) * rhos + ( v_ip(ji,jj,:) + v_il(ji,jj,:) ) * rhow )
pdiag_s(ji,jj) = SUM( sv_i(ji,jj,:) * rhoi )
pdiag_t(ji,jj) = SUM( SUM( e_i(ji,jj,:,:), dim=2 ) ) + SUM( SUM( e_s(ji,jj,:,:), dim=2 ) )
! mass flux
pdiag_fv = wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + &
& wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr
pdiag_fv(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) &
& + wfx_opw(ji,jj) + wfx_res(ji,jj) + wfx_dyn(ji,jj) + wfx_lam(ji,jj) + wfx_pnd (ji,jj) &
& + wfx_snw_sni(ji,jj) + wfx_snw_sum(ji,jj) + wfx_snw_dyn(ji,jj) &
& + wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj) + wfx_spr(ji,jj)
! salt flux
pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam
pdiag_fs(ji,jj) = sfx_bri(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_sub(ji,jj) + sfx_lam(ji,jj)
! heat flux
pdiag_ft = hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr
pdiag_ft(ji,jj) = hfx_sum(ji,jj) + hfx_bom(ji,jj) + hfx_bog(ji,jj) + hfx_dif(ji,jj) + hfx_opw(ji,jj) + hfx_snw(ji,jj) &
& - hfx_thd(ji,jj) - hfx_dyn(ji,jj) - hfx_res(ji,jj) - hfx_sub(ji,jj) - hfx_spr(ji,jj)
END_2D
ELSEIF( icount == 1 ) THEN
! -- mass diag -- !
zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) - pdiag_v ) * r1_Dt_ice &
& + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + &
& wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) &
& - pdiag_fv
DO_2D( 0, 0, 0, 0 )
zdiag_mass(ji,jj) = ( SUM( v_i(ji,jj,:) * rhoi + v_s(ji,jj,:) * rhos &
& + ( v_ip(ji,jj,:) + v_il(ji,jj,:) ) * rhow ) - pdiag_v(ji,jj) ) * r1_Dt_ice &
& + ( wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) + wfx_opw(ji,jj) &
& + wfx_res(ji,jj) + wfx_dyn(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) &
& + wfx_snw_sni(ji,jj) + wfx_snw_sum(ji,jj) + wfx_snw_dyn(ji,jj) &
& + wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj) + wfx_spr(ji,jj) ) &
& - pdiag_fv(ji,jj)
END_2D
IF( MAXVAL( ABS(zdiag_mass) ) > rchk_m * rn_icechk_cel ) ll_stop_m = .TRUE.
!
! -- salt diag -- !
zdiag_salt = ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_Dt_ice &
& + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) &
& - pdiag_fs
DO_2D( 0, 0, 0, 0 )
zdiag_salt(ji,jj) = ( SUM( sv_i(ji,jj,:) * rhoi ) - pdiag_s(ji,jj) ) * r1_Dt_ice &
& + ( sfx_bri(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_sub(ji,jj) + sfx_lam(ji,jj) ) &
& - pdiag_fs(ji,jj)
END_2D
IF( MAXVAL( ABS(zdiag_salt) ) > rchk_s * rn_icechk_cel ) ll_stop_s = .TRUE.
!
! -- heat diag -- !
zdiag_heat = ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) - pdiag_t ) * r1_Dt_ice &
& + ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) &
& - pdiag_ft
DO_2D( 0, 0, 0, 0 )
zdiag_heat(ji,jj) = ( SUM( SUM( e_i(ji,jj,:,:), dim=2 ) ) &
& + SUM( SUM( e_s(ji,jj,:,:), dim=2 ) ) - pdiag_t(ji,jj) ) * r1_Dt_ice &
& + ( hfx_sum(ji,jj) + hfx_bom(ji,jj) + hfx_bog(ji,jj) &
& + hfx_dif(ji,jj) + hfx_opw(ji,jj) + hfx_snw(ji,jj) &
& - hfx_thd(ji,jj) - hfx_dyn(ji,jj) - hfx_res(ji,jj) - hfx_sub(ji,jj) - hfx_spr(ji,jj) ) &
& - pdiag_ft(ji,jj)
END_2D
IF( MAXVAL( ABS(zdiag_heat) ) > rchk_t * rn_icechk_cel ) ll_stop_t = .TRUE.
!
! -- other diags -- !
! a_i < 0
zdiag_amin(:,:) = 0._wp
DO jl = 1, jpl
WHERE( a_i(:,:,jl) < 0._wp ) zdiag_amin(:,:) = 1._wp
WHERE( a_i(A2D(0),jl) < 0._wp ) zdiag_amin(:,:) = 1._wp
ENDDO
! v_i < 0
zdiag_vmin(:,:) = 0._wp
DO jl = 1, jpl
WHERE( v_i(:,:,jl) < 0._wp ) zdiag_vmin(:,:) = 1._wp
WHERE( v_i(A2D(0),jl) < 0._wp ) zdiag_vmin(:,:) = 1._wp
ENDDO
! s_i < 0
zdiag_smin(:,:) = 0._wp
DO jl = 1, jpl
WHERE( s_i(:,:,jl) < 0._wp ) zdiag_smin(:,:) = 1._wp
WHERE( s_i(A2D(0),jl) < 0._wp ) zdiag_smin(:,:) = 1._wp
ENDDO
! e_i < 0
zdiag_emin(:,:) = 0._wp
DO jl = 1, jpl
DO jk = 1, nlay_i
WHERE( e_i(:,:,jk,jl) < 0._wp ) zdiag_emin(:,:) = 1._wp
WHERE( e_i(A2D(0),jk,jl) < 0._wp ) zdiag_emin(:,:) = 1._wp
ENDDO
ENDDO
! a_i > amax
......@@ -528,8 +565,8 @@ CONTAINS
INTEGER :: jl, ji, jj
!!-------------------------------------------------------------------
DO ji = mi0(ki), mi1(ki)
DO jj = mj0(kj), mj1(kj)
DO ji = mi0(ki,nn_hls), mi1(ki,nn_hls)
DO jj = mj0(kj,nn_hls), mj1(kj,nn_hls)
WRITE(numout,*) ' time step ',kt,' ',cd1 ! print title
......@@ -733,10 +770,10 @@ CONTAINS
CALL prt_ctl_info(' ')
CALL prt_ctl_info(' - Stresses : ')
CALL prt_ctl_info(' ~~~~~~~~~~ ')
CALL prt_ctl(tab2d_1=utau , clinfo1= ' utau : ', mask1 = umask, &
& tab2d_2=vtau , clinfo2= ' vtau : ', mask2 = vmask)
CALL prt_ctl(tab2d_1=utau_ice , clinfo1= ' utau_ice : ', mask1 = umask, &
& tab2d_2=vtau_ice , clinfo2= ' vtau_ice : ', mask2 = vmask)
CALL prt_ctl(tab2d_1=utau , clinfo1= ' utau : ', mask1 = tmask, &
& tab2d_2=vtau , clinfo2= ' vtau : ', mask2 = tmask)
CALL prt_ctl(tab2d_1=utau_ice , clinfo1= ' utau_ice : ', mask1 = tmask, &
& tab2d_2=vtau_ice , clinfo2= ' vtau_ice : ', mask2 = tmask)
END SUBROUTINE ice_prt3D
......@@ -751,8 +788,9 @@ CONTAINS
!!-------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ice time-step index
!
REAL(wp), DIMENSION(jpi,jpj,6) :: ztmp
REAL(wp), DIMENSION(6) :: zchk
INTEGER :: ji, jj ! dummy loop index
REAL(wp), DIMENSION(A2D(0),6) :: ztmp
REAL(wp), DIMENSION(6) :: zchk
!!-------------------------------------------------------------------
!
IF( kt == nit000 .AND. lwp ) THEN
......@@ -762,25 +800,27 @@ CONTAINS
ENDIF
!
! -- 2D budgets (must be close to 0) -- !
ztmp(:,:,1) = wfx_ice (:,:) + wfx_snw (:,:) + wfx_spr (:,:) + wfx_sub(:,:) + wfx_pnd(:,:) &
& + diag_vice(:,:) + diag_vsnw(:,:) + diag_vpnd(:,:) - diag_adv_mass(:,:)
ztmp(:,:,2) = sfx(:,:) + diag_sice(:,:) - diag_adv_salt(:,:)
ztmp(:,:,3) = qt_oce_ai(:,:) - qt_atm_oi(:,:) + diag_heat(:,:) - diag_adv_heat(:,:)
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,1) = wfx_ice (ji,jj) + wfx_snw (ji,jj) + wfx_spr (ji,jj) + wfx_sub(ji,jj) + wfx_pnd(ji,jj) &
& + diag_vice(ji,jj) + diag_vsnw(ji,jj) + diag_vpnd(ji,jj) - diag_adv_mass(ji,jj)
ztmp(ji,jj,2) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj)
ztmp(ji,jj,3) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj)
END_2D
! write outputs
CALL iom_put( 'icedrift_mass', ztmp(:,:,1) )
CALL iom_put( 'icedrift_salt', ztmp(:,:,2) )
CALL iom_put( 'icedrift_heat', ztmp(:,:,3) )
! -- 1D budgets -- !
ztmp(:,:,1) = ztmp(:,:,1) * e1e2t * rDt_ice ! mass
ztmp(:,:,2) = ztmp(:,:,2) * e1e2t * rDt_ice * 1.e-3 ! salt
ztmp(:,:,3) = ztmp(:,:,3) * e1e2t ! heat
ztmp(:,:,4) = diag_adv_mass * e1e2t * rDt_ice
ztmp(:,:,5) = diag_adv_salt * e1e2t * rDt_ice * 1.e-3
ztmp(:,:,6) = diag_adv_heat * e1e2t
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,1) = ztmp(ji,jj,1) * e1e2t(ji,jj) * rDt_ice ! mass
ztmp(ji,jj,2) = ztmp(ji,jj,2) * e1e2t(ji,jj) * rDt_ice * 1.e-3 ! salt
ztmp(ji,jj,3) = ztmp(ji,jj,3) * e1e2t(ji,jj) ! heat
ztmp(ji,jj,4) = diag_adv_mass(ji,jj) * e1e2t(ji,jj) * rDt_ice
ztmp(ji,jj,5) = diag_adv_salt(ji,jj) * e1e2t(ji,jj) * rDt_ice * 1.e-3
ztmp(ji,jj,6) = diag_adv_heat(ji,jj) * e1e2t(ji,jj)
END_2D
! global sums
zchk(1:6) = glob_sum_vec( 'icectl', ztmp(:,:,1:6) )
......
src/ICE/icedia.F90
View file @ 3efbcb17
......@@ -33,6 +33,9 @@ MODULE icedia
PUBLIC ice_dia ! called by icestp.F90
PUBLIC ice_dia_init ! called in icestp.F90
!! * Substitutions
# include "do_loop_substitute.h90"
REAL(wp), SAVE :: r1_area ! inverse of the ocean area
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents
REAL(wp) :: frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot ! global forcing trends
......@@ -48,7 +51,7 @@ CONTAINS
!!---------------------------------------------------------------------!
!! *** ROUTINE ice_dia_alloc ***
!!---------------------------------------------------------------------!
ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ice_dia_alloc )
ALLOCATE( vol_loc_ini(A2D(0)), sal_loc_ini(A2D(0)), tem_loc_ini(A2D(0)), STAT=ice_dia_alloc )
CALL mpp_sum ( 'icedia', ice_dia_alloc )
IF( ice_dia_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_dia_alloc: failed to allocate arrays' )
......@@ -64,8 +67,9 @@ CONTAINS
!!---------------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time step
!!
REAL(wp), DIMENSION(jpi,jpj,16) :: ztmp
REAL(wp), DIMENSION(16) :: zbg
INTEGER :: ji, jj ! dummy loop index
REAL(wp), DIMENSION(A2D(0),16) :: ztmp
REAL(wp), DIMENSION(16) :: zbg
!!---------------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('ice_dia')
......@@ -85,31 +89,32 @@ CONTAINS
! 1 - Trends due to forcing !
! ---------------------------!
! they must be kept outside an IF(iom_use) because of the call to dia_rst below
ztmp(:,:,1) = - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ! freshwater flux ice/snow-ocean
ztmp(:,:,2) = - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ! freshwater flux ice/snow-atm
ztmp(:,:,3) = - sfx (:,:) * e1e2t(:,:) ! salt fluxes ice/snow-ocean
ztmp(:,:,4) = qt_atm_oi(:,:) * e1e2t(:,:) ! heat on top of ice-ocean
ztmp(:,:,5) = qt_oce_ai(:,:) * e1e2t(:,:) ! heat on top of ocean (and below ice)
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,1) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) ) * e1e2t(ji,jj) ! freshwater flux ice/snow-ocean
ztmp(ji,jj,2) = - ( wfx_sub(ji,jj) + wfx_spr(ji,jj) ) * e1e2t(ji,jj) ! freshwater flux ice/snow-atm
ztmp(ji,jj,3) = - sfx (ji,jj) * e1e2t(ji,jj) ! salt fluxes ice/snow-ocean
ztmp(ji,jj,4) = qt_atm_oi(ji,jj) * e1e2t(ji,jj) ! heat on top of ice-ocean
ztmp(ji,jj,5) = qt_oce_ai(ji,jj) * e1e2t(ji,jj) ! heat on top of ocean (and below ice)
END_2D
! ----------------------- !
! 2 - Contents !
! ----------------------- !
IF( iom_use('ibgvol_tot' ) ) ztmp(:,:,6 ) = vt_i (:,:) * e1e2t(:,:) ! ice volume
IF( iom_use('sbgvol_tot' ) ) ztmp(:,:,7 ) = vt_s (:,:) * e1e2t(:,:) ! snow volume
IF( iom_use('ibgarea_tot') ) ztmp(:,:,8 ) = at_i (:,:) * e1e2t(:,:) ! area
IF( iom_use('ibgsalt_tot') ) ztmp(:,:,9 ) = st_i (:,:) * e1e2t(:,:) ! salt content
IF( iom_use('ibgheat_tot') ) ztmp(:,:,10) = et_i (:,:) * e1e2t(:,:) ! heat content
IF( iom_use('sbgheat_tot') ) ztmp(:,:,11) = et_s (:,:) * e1e2t(:,:) ! heat content
IF( iom_use('ipbgvol_tot') ) ztmp(:,:,12) = vt_ip(:,:) * e1e2t(:,:) ! ice pond volume
IF( iom_use('ilbgvol_tot') ) ztmp(:,:,13) = vt_il(:,:) * e1e2t(:,:) ! ice pond lid volume
IF( iom_use('ibgvol_tot' ) ) ztmp(:,:,6 ) = vt_i (A2D(0)) * e1e2t(A2D(0)) ! ice volume
IF( iom_use('sbgvol_tot' ) ) ztmp(:,:,7 ) = vt_s (A2D(0)) * e1e2t(A2D(0)) ! snow volume
IF( iom_use('ibgarea_tot') ) ztmp(:,:,8 ) = at_i (A2D(0)) * e1e2t(A2D(0)) ! area
IF( iom_use('ibgsalt_tot') ) ztmp(:,:,9 ) = st_i (:,:) * e1e2t(A2D(0)) ! salt content
IF( iom_use('ibgheat_tot') ) ztmp(:,:,10) = et_i (:,:) * e1e2t(A2D(0)) ! heat content
IF( iom_use('sbgheat_tot') ) ztmp(:,:,11) = et_s (:,:) * e1e2t(A2D(0)) ! heat content
IF( iom_use('ipbgvol_tot') ) ztmp(:,:,12) = vt_ip(A2D(0)) * e1e2t(A2D(0)) ! ice pond volume
IF( iom_use('ilbgvol_tot') ) ztmp(:,:,13) = vt_il(A2D(0)) * e1e2t(A2D(0)) ! ice pond lid volume
! ---------------------------------- !
! 3 - Content variations and drifts !
! ---------------------------------- !
IF( iom_use('ibgvolume') ) ztmp(:,:,14) = ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ! freshwater trend
IF( iom_use('ibgsaltco') ) ztmp(:,:,15) = ( rhoi*st_i(:,:) - sal_loc_ini(:,:) ) * e1e2t(:,:) ! salt content trend
IF( iom_use('ibgvolume') ) ztmp(:,:,14) = ( rhoi*vt_i(A2D(0)) + rhos*vt_s(A2D(0)) - vol_loc_ini(:,:) ) * e1e2t(A2D(0)) ! freshwater trend
IF( iom_use('ibgsaltco') ) ztmp(:,:,15) = ( rhoi*st_i(:,:) - sal_loc_ini(:,:) ) * e1e2t(A2D(0)) ! salt content trend
IF( iom_use('ibgheatco') .OR. iom_use('ibgheatfx') ) &
& ztmp(:,:,16) = ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:) ) * e1e2t(:,:) ! heat content trend
& ztmp(:,:,16) = ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:) ) * e1e2t(A2D(0)) ! heat content trend
! global sum
zbg(1:16) = glob_sum_vec( 'icedia', ztmp(:,:,1:16) )
......@@ -261,9 +266,9 @@ CONTAINS
frc_tembot = 0._wp
frc_sal = 0._wp
! record initial ice volume, salt and temp
vol_loc_ini(:,:) = rhoi * vt_i(:,:) + rhos * vt_s(:,:) ! ice/snow volume (kg/m2)
tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:) ! ice/snow heat content (J)
sal_loc_ini(:,:) = rhoi * st_i(:,:) ! ice salt content (pss*kg/m2)
vol_loc_ini(:,:) = rhoi * vt_i(A2D(0)) + rhos * vt_s(A2D(0)) ! ice/snow volume (kg/m2)
tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:) ! ice/snow heat content (J)
sal_loc_ini(:,:) = rhoi * st_i(:,:) ! ice salt content (pss*kg/m2)
ENDIF
!
ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
......
src/ICE/icedyn.F90
View file @ 3efbcb17
......@@ -93,8 +93,8 @@ CONTAINS
!
! retrieve thickness from volume for landfast param. and UMx advection scheme
WHERE( a_i(:,:,:) >= epsi20 )
h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:)
h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:)
h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:)
h_s(:,:,:) = v_s(:,:,:) / a_i(:,:,:)
ELSEWHERE
h_i(:,:,:) = 0._wp
h_s(:,:,:) = 0._wp
......@@ -162,15 +162,12 @@ CONTAINS
CASE ( np_dynADV1D , np_dynADV2D )
ALLOCATE( zdivu_i(jpi,jpj) )
ALLOCATE( zdivu_i(A2D(0)) )
DO_2D( 0, 0, 0, 0 )
zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
END_2D
CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1.0_wp )
! output
CALL iom_put( 'icediv' , zdivu_i )
DEALLOCATE( zdivu_i )
END SELECT
......
src/ICE/icedyn_adv.F90
View file @ 3efbcb17
......@@ -41,6 +41,8 @@ MODULE icedyn_adv
! ** namelist (namdyn_adv) **
INTEGER :: nn_UMx ! order of the UMx advection scheme
!
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: icedyn_adv.F90 13472 2020-09-16 13:05:19Z smasson $
......@@ -92,11 +94,11 @@ CONTAINS
!------------
! diagnostics
!------------
diag_trp_ei(:,:) = SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice
diag_trp_es(:,:) = SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice
diag_trp_sv(:,:) = SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice
diag_trp_vi(:,:) = SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice
diag_trp_vs(:,:) = SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice
diag_trp_ei(:,:) = SUM(SUM( e_i (A2D(0),1:nlay_i,:) - e_i_b (A2D(0),1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice
diag_trp_es(:,:) = SUM(SUM( e_s (A2D(0),1:nlay_s,:) - e_s_b (A2D(0),1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice
diag_trp_sv(:,:) = SUM( sv_i(A2D(0),:) - sv_i_b(A2D(0),:) , dim=3 ) * r1_Dt_ice
diag_trp_vi(:,:) = SUM( v_i (A2D(0),:) - v_i_b (A2D(0),:) , dim=3 ) * r1_Dt_ice
diag_trp_vs(:,:) = SUM( v_s (A2D(0),:) - v_s_b (A2D(0),:) , dim=3 ) * r1_Dt_ice
IF( iom_use('icemtrp') ) CALL iom_put( 'icemtrp' , diag_trp_vi * rhoi ) ! ice mass transport
IF( iom_use('snwmtrp') ) CALL iom_put( 'snwmtrp' , diag_trp_vs * rhos ) ! snw mass transport
IF( iom_use('salmtrp') ) CALL iom_put( 'salmtrp' , diag_trp_sv * rhoi * 1.e-03 ) ! salt mass transport (kg/m2/s)
......
src/ICE/icedyn_adv_pra.F90
View file @ 3efbcb17
......@@ -89,7 +89,7 @@ CONTAINS
INTEGER :: icycle ! number of sub-timestep for the advection
REAL(wp) :: zdt, z1_dt ! - -
REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication
REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2
REAL(wp), DIMENSION(A2D(0)) :: zati1, zati2
REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max, zs_i, zsi_max
REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: ze_i, zei_max
......@@ -100,7 +100,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es
REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei
!! diagnostics
REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
REAL(wp), DIMENSION(A2D(0)) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
!!----------------------------------------------------------------------
!
IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
......@@ -129,8 +129,7 @@ CONTAINS
END DO
CALL icemax4D( ze_i , zei_max )
CALL icemax4D( ze_s , zes_max )
CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp, zes_max, 'T', 1._wp )
!
!
! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
......@@ -155,14 +154,14 @@ CONTAINS
DO jt = 1, icycle
! diagnostics
zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow
zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi
zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 )
zdiag_adv_mass(:,:) = SUM( pv_i (A2D(0),:) , dim=3 ) * rhoi + SUM( pv_s (A2D(0),:) , dim=3 ) * rhos &
& + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow
zdiag_adv_salt(:,:) = SUM( psv_i(A2D(0),:) , dim=3 ) * rhoi
zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(A2D(0),1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(A2D(0),1:nlay_s,:) , dim=4 ), dim=3 )
! record at_i before advection (for open water)
zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
zati1(:,:) = SUM( pa_i(A2D(0),:), dim=3 )
! --- transported fields --- !
DO jl = 1, jpl
......@@ -275,8 +274,8 @@ CONTAINS
& , z0oi , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp & ! ice age
& , sxxage, 'T', 1._wp, syyage, 'T', 1._wp, sxyage, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_pra', z0es , 'T', 1._wp, sxc0 , 'T', -1._wp, syc0 , 'T', -1._wp & ! snw enthalpy
& , sxxc0 , 'T', 1._wp, syyc0 , 'T', 1._wp, sxyc0 , 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_pra', z0ei , 'T', 1._wp, sxe , 'T', -1._wp, sye , 'T', -1._wp & ! ice enthalpy
& , sxxc0 , 'T', 1._wp, syyc0 , 'T', 1._wp, sxyc0 , 'T', 1._wp &
& , z0ei , 'T', 1._wp, sxe , 'T', -1._wp, sye , 'T', -1._wp & ! ice enthalpy
& , sxxe , 'T', 1._wp, syye , 'T', 1._wp, sxye , 'T', 1._wp )
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
IF( ln_pnd_lids ) THEN
......@@ -317,7 +316,7 @@ CONTAINS
END DO
!
! derive open water from ice concentration
zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
zati2(:,:) = SUM( pa_i(A2D(0),:), dim=3 )
DO_2D( 0, 0, 0, 0 )
pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water
& - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
......@@ -325,13 +324,13 @@ CONTAINS
CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1.0_wp )
!
! --- diagnostics --- !
diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow &
diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (A2D(0),:) , dim=3 ) * rhoi + SUM( pv_s (A2D(0),:) , dim=3 ) * rhos &
& + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow &
& - zdiag_adv_mass(:,:) ) * z1_dt
diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi &
diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(A2D(0),:) , dim=3 ) * rhoi &
& - zdiag_adv_salt(:,:) ) * z1_dt
diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) &
diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(A2D(0),1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(A2D(0),1:nlay_s,:) , dim=4 ), dim=3 ) &
& - zdiag_adv_heat(:,:) ) * z1_dt
!
! --- Ensure non-negative fields --- !
......
src/ICE/icedyn_adv_umx.F90
View file @ 3efbcb17
......@@ -104,7 +104,7 @@ CONTAINS
!
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs
!! diagnostics
REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
REAL(wp), DIMENSION(A2D(0)) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
!!----------------------------------------------------------------------
!
IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
......@@ -133,8 +133,7 @@ CONTAINS
END DO
CALL icemax4D( ze_i , zei_max )
CALL icemax4D( ze_s , zes_max )
CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp, zes_max, 'T', 1._wp )
!
!
! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
......@@ -182,11 +181,11 @@ CONTAINS
DO jt = 1, icycle
! diagnostics
zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow
zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi
zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 )
zdiag_adv_mass(:,:) = SUM( pv_i (A2D(0),:) , dim=3 ) * rhoi + SUM( pv_s (A2D(0),:) , dim=3 ) * rhos &
& + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow
zdiag_adv_salt(:,:) = SUM( psv_i(A2D(0),:) , dim=3 ) * rhoi
zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(A2D(0),1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(A2D(0),1:nlay_s,:) , dim=4 ), dim=3 )
! record at_i before advection (for open water)
zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
......@@ -370,8 +369,7 @@ CONTAINS
ELSE
CALL lbc_lnk( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp )
ENDIF
CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp, pe_s, 'T', 1._wp )
!
!== Open water area ==!
zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
......@@ -382,13 +380,13 @@ CONTAINS
CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp )
!
! --- diagnostics --- !
diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow &
diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (A2D(0),:) , dim=3 ) * rhoi + SUM( pv_s (A2D(0),:) , dim=3 ) * rhos &
& + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow &
& - zdiag_adv_mass(:,:) ) * z1_dt
diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi &
diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(A2D(0),:) , dim=3 ) * rhoi &
& - zdiag_adv_salt(:,:) ) * z1_dt
diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) &
diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(A2D(0),1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(A2D(0),1:nlay_s,:) , dim=4 ), dim=3 ) &
& - zdiag_adv_heat(:,:) ) * z1_dt
!
! --- Ensure non-negative fields and in-bound thicknesses --- !
......
src/ICE/icedyn_rdgrft.F90
View file @ 3efbcb17
......@@ -174,7 +174,7 @@ CONTAINS
!
npti = 0 ; nptidx(:) = 0
ipti = 0 ; iptidx(:) = 0
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF ( at_i(ji,jj) > epsi10 ) THEN
npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji
......@@ -280,8 +280,15 @@ CONTAINS
CALL ice_dyn_1d2d( 2 ) ! --- Move to 2D arrays --- !
ENDIF
! clem: the 3 lbc below could be avoided if calculations above were performed over the full domain
! but we think it is more efficient this way => to check?
CALL lbc_lnk( 'icedyn_rdgrft', ato_i , 'T', 1._wp )
CALL lbc_lnk( 'icedyn_rdgrft', a_i , 'T', 1._wp, v_i , 'T', 1._wp, v_s , 'T', 1._wp, sv_i, 'T', 1._wp, oa_i, 'T', 1._wp, &
& a_ip , 'T', 1._wp, v_ip, 'T', 1._wp, v_il, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_rdgrft', e_i , 'T', 1._wp, e_s , 'T', 1._wp )
CALL ice_var_agg( 1 )
! clem: I think we can comment this line but I am not sure it does not change results
!!$ CALL ice_var_agg( 1 )
! controls
IF( sn_cfctl%l_prtctl ) CALL ice_prt3D('icedyn_rdgrft') ! prints
......@@ -302,8 +309,9 @@ CONTAINS
!! ** Method : Compute the thickness distribution of the ice and open water
!! participating in ridging and of the resulting ridges.
!!-------------------------------------------------------------------
REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i, pclosing_net
REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i
REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i
REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i
REAL(wp), DIMENSION(:) , INTENT(in), OPTIONAL :: pclosing_net
!!
INTEGER :: ji, jl ! dummy loop indices
REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar
......@@ -504,39 +512,43 @@ CONTAINS
END DO
END DO
!
! 3) closing_gross
!-----------------
! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.
! NOTE: 0 < aksum <= 1
WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
ELSEWHERE ; closing_gross(1:npti) = 0._wp
END WHERE
! correction to closing rate if excessive ice removal
!----------------------------------------------------
! Reduce the closing rate if more than 100% of any ice category would be removed
! Reduce the opening rate in proportion
DO jl = 1, jpl
IF( PRESENT( pclosing_net ) ) THEN
!
! 3) closing_gross
!-----------------
! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.
! NOTE: 0 < aksum <= 1
WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
ELSEWHERE ; closing_gross(1:npti) = 0._wp
END WHERE
! correction to closing rate if excessive ice removal
!----------------------------------------------------
! Reduce the closing rate if more than 100% of any ice category would be removed
! Reduce the opening rate in proportion
DO jl = 1, jpl
DO ji = 1, npti
zfac = apartf(ji,jl) * closing_gross(ji) * rDt_ice
IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice
ENDIF
END DO
END DO
! 4) correction to opening if excessive open water removal
!---------------------------------------------------------
! Reduce the closing rate if more than 100% of the open water would be removed
! Reduce the opening rate in proportion
DO ji = 1, npti
zfac = apartf(ji,jl) * closing_gross(ji) * rDt_ice
IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice
zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice
IF( zfac < 0._wp ) THEN ! would lead to negative ato_i
opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice
ELSEIF( zfac > zasum(ji) ) THEN ! would lead to ato_i > asum
opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice
ENDIF
END DO
END DO
! 4) correction to opening if excessive open water removal
!---------------------------------------------------------
! Reduce the closing rate if more than 100% of the open water would be removed
! Reduce the opening rate in proportion
DO ji = 1, npti
zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice
IF( zfac < 0._wp ) THEN ! would lead to negative ato_i
opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice
ELSEIF( zfac > zasum(ji) ) THEN ! would lead to ato_i > asum
opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice
ENDIF
END DO
!
ENDIF
!
END SUBROUTINE rdgrft_prep
......@@ -861,7 +873,8 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jl ! dummy loop indices
REAL(wp) :: z1_3 ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: zmsk, zworka ! temporary array used here
REAL(wp), DIMENSION(A2D(0)) :: zworka ! temporary array used here
REAL(wp), DIMENSION(jpi,jpj) :: zmsk ! temporary array used here
!!
LOGICAL :: ln_str_R75
REAL(wp) :: zhi, zcp
......@@ -886,7 +899,7 @@ CONTAINS
!
! Identify grid cells with ice
npti = 0 ; nptidx(:) = 0
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF ( at_i(ji,jj) > epsi10 ) THEN
npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji
......@@ -899,7 +912,7 @@ CONTAINS
CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i )
CALL tab_2d_1d( npti, nptidx(1:npti), zstrength(1:npti) , strength )
CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d )
!
zaksum(1:npti) = apartf(1:npti,0) !clem: aksum should be defined in the header => local to module
DO jl = 1, jpl
......@@ -956,12 +969,20 @@ CONTAINS
CALL tab_1d_2d( npti, nptidx(1:npti), zstrength(1:npti), strength )
!
ENDIF
CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp ) ! this call could be removed if calculations were done on the full domain
! ! but we decided it is more efficient this way
!
CASE ( np_strh79 ) !== Hibler(1979)'s method ==!
strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - at_i(:,:) ) ) * zmsk(:,:)
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
strength(ji,jj) = rn_pstar * SUM( v_i(ji,jj,:) ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) * zmsk(ji,jj)
END_2D
!
CASE ( np_strcst ) !== Constant strength ==!
strength(:,:) = rn_str * zmsk(:,:)
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
strength(ji,jj) = rn_str * zmsk(ji,jj)
END_2D
!
END SELECT
!
......
src/ICE/icedyn_rhg_eap.F90
View file @ 3efbcb17
......@@ -125,12 +125,12 @@ CONTAINS
!! Bouillon et al., Ocean Modelling 2013
!! Kimmritz et al., Ocean Modelling 2016 & 2017
!!-------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i !
REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i !
REAL(wp), DIMENSION(:,:), INTENT(inout) :: paniso_11 , paniso_12 ! structure tensor components
REAL(wp), DIMENSION(:,:), INTENT(inout) :: prdg_conv ! for ridging
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i !
REAL(wp), DIMENSION(A2D(0)), INTENT( out) :: pshear_i , pdivu_i , pdelta_i !
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: paniso_11 , paniso_12 ! structure tensor components
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: prdg_conv ! for ridging
!!
INTEGER :: ji, jj ! dummy loop indices
INTEGER :: jter ! local integers
......@@ -148,15 +148,15 @@ CONTAINS
!
REAL(wp) :: zintb, zintn ! dummy argument
REAL(wp) :: zfac_x, zfac_y
REAL(wp) :: zshear, zdum1, zdum2
REAL(wp) :: zdum1, zdum2
REAL(wp) :: zstressptmp, zstressmtmp, zstress12tmpF ! anisotropic stress tensor components
REAL(wp) :: zalphar, zalphas ! for mechanical redistribution
REAL(wp) :: zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A ! for structure tensor evolution
!
REAL(wp), DIMENSION(jpi,jpj) :: zstress12tmp ! anisotropic stress tensor component for regridding
REAL(wp), DIMENSION(jpi,jpj) :: zyield11, zyield22, zyield12 ! yield surface tensor for history
REAL(wp), DIMENSION(A2D(0)) :: zyield11, zyield22, zyield12 ! yield surface tensor for history
REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points
REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension
REAL(wp), DIMENSION(A2D(0)) :: zten_i, zshear ! tension, shear
REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017
!
REAL(wp), DIMENSION(jpi,jpj) :: zdt_m ! (dt / ice-snow_mass) on T points
......@@ -178,7 +178,6 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast)
REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast)
!
REAL(wp), DIMENSION(jpi,jpj) :: zmsk00, zmsk15
REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays
REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence
......@@ -186,7 +185,8 @@ CONTAINS
REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small
REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small
!! --- check convergence
REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice
REAL(wp), DIMENSION(A2D(0)) :: zmsk00, zmsk15
REAL(wp), DIMENSION(A2D(0)) :: zu_ice, zv_ice
!! --- diags
REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p
......@@ -202,11 +202,11 @@ CONTAINS
IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology'
!
! for diagnostics and convergence tests
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice
END_2D
IF( nn_rhg_chkcvg > 0 ) THEN
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
END_2D
ENDIF
......@@ -283,7 +283,14 @@ CONTAINS
! non-embedded sea ice: use ocean surface for slope calculation
zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
DO_2D( 0, 0, 0, 0 )
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zm1 = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) ) ! Ice/snow mass at U-V points
!!$ zm1 = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) + rhow * (vt_ip(ji,jj) + vt_il(ji,jj)) ) ! clem: this should replace the above
zmf (ji,jj) = zm1 * ff_t(ji,jj) ! Coriolis at T points (m*f)
zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) ! dt/m at T points (for alpha and beta coefficients)
END_2D
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
! ice fraction at U-V points
zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
......@@ -291,8 +298,11 @@ CONTAINS
! Ice/snow mass at U-V points
zm1 = ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ) )
!!$ zm1 = ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ) + rhow * (vt_ip(ji ,jj ) + vt_il(ji ,jj )) ) ! clem: this should replace the above
zm2 = ( rhos * vt_s(ji+1,jj ) + rhoi * vt_i(ji+1,jj ) )
!!$ zm2 = ( rhos * vt_s(ji+1,jj ) + rhoi * vt_i(ji+1,jj ) + rhow * (vt_ip(ji+1,jj ) + vt_il(ji+1,jj )) ) ! clem: this should replace the above
zm3 = ( rhos * vt_s(ji ,jj+1) + rhoi * vt_i(ji ,jj+1) )
!!$ zm3 = ( rhos * vt_s(ji ,jj+1) + rhoi * vt_i(ji ,jj+1) + rhow * (vt_ip(ji ,jj+1) + vt_il(ji ,jj+1)) ) ! clem: this should replace the above
zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
......@@ -300,19 +310,17 @@ CONTAINS
v_oceU(ji,jj) = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
u_oceV(ji,jj) = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
! Coriolis at T points (m*f)
zmf(ji,jj) = zm1 * ff_t(ji,jj)
! dt/m at T points (for alpha and beta coefficients)
zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin )
! m/dt
zmU_t(ji,jj) = zmassU * z1_dtevp
zmV_t(ji,jj) = zmassV * z1_dtevp
! Drag ice-atm.
ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
! and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
ztaux_ai(ji,jj) = zaU(ji,jj) * 0.5_wp * ( utau_ice(ji,jj) + utau_ice(ji+1,jj) ) * &
& ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
ztauy_ai(ji,jj) = zaV(ji,jj) * 0.5_wp * ( vtau_ice(ji,jj) + vtau_ice(ji,jj+1) ) * &
& ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
......@@ -329,12 +337,11 @@ CONTAINS
ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp )
!
! !== Landfast ice parameterization ==!
!
IF( ln_landfast_L16 ) THEN !-- Lemieux 2016
DO_2D( 0, 0, 0, 0 )
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
! ice thickness at U-V points
zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
......@@ -348,10 +355,9 @@ CONTAINS
zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', tau_icebfr(:,:), 'T', 1.0_wp )
!
ELSE !-- no landfast
DO_2D( 0, 0, 0, 0 )
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
ztaux_base(ji,jj) = 0._wp
ztauy_base(ji,jj) = 0._wp
END_2D
......@@ -368,14 +374,14 @@ CONTAINS
!
! convergence test
IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
END_2D
ENDIF
! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
DO_2D( 1, 0, 1, 0 )
DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
! shear at F points
zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &
......@@ -406,15 +412,13 @@ CONTAINS
! delta at T points
zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1.0_wp )
! P/delta at T points
DO_2D( 1, 1, 1, 1 )
! P/delta at T points
zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1.0_wp, zp_delt, 'T', 1.0_wp )
DO_2D( 0, 1, 0, 1 ) ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
DO_2D( 0, 0, 0, 0 )
! shear at T points
zdsT = ( zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
......@@ -467,16 +471,17 @@ CONTAINS
zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zstressptmp ) * z1_alph1
zs2(ji,jj) = ( zs2(ji,jj) * zalph1 + zstressmtmp ) * z1_alph1
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', zstress12tmp, 'T', 1.0_wp , paniso_11, 'T', 1.0_wp , paniso_12, 'T', 1.0_wp)
CALL lbc_lnk( 'icedyn_rhg_eap', zstress12tmp, 'T', 1.0_wp , paniso_11, 'T', 1.0_wp , paniso_12, 'T', 1.0_wp, &
& zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp )
! Save beta at T-points for further computations
IF( ln_aEVP ) THEN
DO_2D( 1, 1, 1, 1 )
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
END_2D
ENDIF
DO_2D( 1, 0, 1, 0 )
DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
! stress12tmp at F points
zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1) &
& + zstress12tmp(ji,jj ) * e1e2t(ji,jj ) + zstress12tmp(ji+1,jj ) * e1e2t(ji+1,jj ) &
......@@ -495,10 +500,9 @@ CONTAINS
zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )
! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
DO_2D( 0, 0, 0, 0 )
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
! !--- U points
zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
& + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &
......@@ -526,7 +530,7 @@ CONTAINS
! Bouillon et al. 2009 (eq 34-35) => stable
IF( MOD(jter,2) == 0 ) THEN ! even iterations
!
DO_2D( 0, 0, 0, 0 )
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
! !--- tau_io/(v_oce - v_ice)
zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
......@@ -570,13 +574,7 @@ CONTAINS
& ) * zmsk00y(ji,jj)
ENDIF
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp )
!
#if defined key_agrif
!! CALL agrif_interp_ice( 'V', jter, nn_nevp )
CALL agrif_interp_ice( 'V' )
#endif
IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp )
!
DO_2D( 0, 0, 0, 0 )
! !--- tau_io/(u_oce - u_ice)
......@@ -622,17 +620,13 @@ CONTAINS
& ) * zmsk00x(ji,jj)
ENDIF
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp )
!
#if defined key_agrif
!! CALL agrif_interp_ice( 'U', jter, nn_nevp )
CALL agrif_interp_ice( 'U' )
#endif
IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
IF( nn_hls == 1 ) THEN ; CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp )
ELSE ; CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp )
ENDIF
!
ELSE ! odd iterations
!
DO_2D( 0, 0, 0, 0 )
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
! !--- tau_io/(u_oce - u_ice)
zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
......@@ -676,13 +670,7 @@ CONTAINS
& ) * zmsk00x(ji,jj)
ENDIF
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp )
!
#if defined key_agrif
!! CALL agrif_interp_ice( 'U', jter, nn_nevp )
CALL agrif_interp_ice( 'U' )
#endif
IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp )
!
DO_2D( 0, 0, 0, 0 )
! !--- tau_io/(v_oce - v_ice)
......@@ -728,15 +716,19 @@ CONTAINS
& ) * zmsk00y(ji,jj)
ENDIF
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp )
IF( nn_hls == 1 ) THEN ; CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp )
ELSE ; CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp )
ENDIF
!
ENDIF
#if defined key_agrif
!! CALL agrif_interp_ice( 'V', jter, nn_nevp )
CALL agrif_interp_ice( 'V' )
!! CALL agrif_interp_ice( 'U', jter, nn_nevp )
!! CALL agrif_interp_ice( 'V', jter, nn_nevp )
CALL agrif_interp_ice( 'U' )
CALL agrif_interp_ice( 'V' )
#endif
IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
!
ENDIF
IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
! convergence test
IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 )
......@@ -751,7 +743,7 @@ CONTAINS
!------------------------------------------------------------------------------!
! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
!------------------------------------------------------------------------------!
DO_2D( 1, 0, 1, 0 )
DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
! shear at F points
zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &
......@@ -775,22 +767,30 @@ CONTAINS
& + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &
& ) * 0.25_wp * r1_e1e2t(ji,jj)
! shear at T points
! maximum shear rate at T points (includes tension, output only)
pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
! shear at T-points
zshear(ji,jj) = SQRT( zds2 )
! divergence at T points
pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
& ) * r1_e1e2t(ji,jj)
! delta at T points
zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta
rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
zdelta(ji,jj) = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta
! delta at T points
rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0
pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch
! it seems that deformation used for advection and mech redistribution is delta*
! MV in principle adding creep limit is a regularization for viscosity not for delta
! delta_star should not (in my view) be used in a replacement for delta
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp, &
& zten_i, 'T', 1.0_wp, zs1 , 'T', 1.0_wp, zs2 , 'T', 1.0_wp, &
& zs1, 'T', 1.0_wp, zs2 , 'T', 1.0_wp, &
& zs12, 'F', 1.0_wp )
! --- Store the stress tensor for the next time step --- !
......@@ -803,45 +803,38 @@ CONTAINS
! 5) diagnostics
!------------------------------------------------------------------------------!
! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
& iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
!
CALL lbc_lnk( 'icedyn_rhg_eap', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, ztaux_ai, 'U', -1.0_wp, &
& ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp )
!
CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
ENDIF
IF( iom_use('utau_oi') ) CALL iom_put( 'utau_oi' , ztaux_oi(A2D(0)) * zmsk00 )
IF( iom_use('vtau_oi') ) CALL iom_put( 'vtau_oi' , ztauy_oi(A2D(0)) * zmsk00 )
IF( iom_use('utau_ai') ) CALL iom_put( 'utau_ai' , ztaux_ai(A2D(0)) * zmsk00 )
IF( iom_use('vtau_ai') ) CALL iom_put( 'vtau_ai' , ztauy_ai(A2D(0)) * zmsk00 )
IF( iom_use('utau_bi') ) CALL iom_put( 'utau_bi' , ztaux_bi(A2D(0)) * zmsk00 )
IF( iom_use('vtau_bi') ) CALL iom_put( 'vtau_bi' , ztauy_bi(A2D(0)) * zmsk00 )
! --- divergence, shear and strength --- !
IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence
IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear
IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta
IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength
IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i (A2D(0)) * zmsk00 ) ! divergence
IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i(A2D(0)) * zmsk00 ) ! shear
IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength(A2D(0)) * zmsk00 ) ! strength
IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , zdelta (A2D(0)) * zmsk00 ) ! delta
! --- Stress tensor invariants (SIMIP diags) --- !
IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
!
ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
ALLOCATE( zsig_I(A2D(0)) , zsig_II(A2D(0)) )
!
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
! Ice stresses
! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
! These are NOT stress tensor components, neither stress invariants, neither stress principal components
! I know, this can be confusing...
zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )
zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
zfac = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) ! viscosity
zsig1 = zfac * ( pdivu_i(ji,jj) - zdelta(ji,jj) )
zsig2 = zfac * z1_ecc2 * zten_i(ji,jj)
zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)
zsig12 = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp
! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure
zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress
zsig_I (ji,jj) = 0.5_wp * zsig1
zsig_II(ji,jj) = 0.5_wp * SQRT ( zsig2 * zsig2 + 4._wp * zsig12 * zsig12 )
END_2D
!
......@@ -859,21 +852,20 @@ CONTAINS
! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
!
ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
ALLOCATE( zsig1_p(A2D(0)) , zsig2_p(A2D(0)) , zsig_I(A2D(0)) , zsig_II(A2D(0)) )
!
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
! and **deformations** at current iterates
! For EVP solvers, ice stresses at current iterates can be used
! following Lemieux & Dupont (2020)
zfac = zp_delt(ji,jj)
zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
zfac = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
zsig1 = zfac * ( pdivu_i(ji,jj) - zdelta(ji,jj) )
zsig2 = zfac * z1_ecc2 * zten_i(ji,jj)
zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)
zsig12 = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp
! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure
zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress
zsig_I(ji,jj) = 0.5_wp * zsig1 ! normal stress
zsig_II(ji,jj) = 0.5_wp * SQRT ( zsig2 * zsig2 + 4._wp * zsig12 * zsig12 ) ! max shear stress
! Normalized principal stresses (used to display the ellipse)
z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) )
......@@ -881,8 +873,8 @@ CONTAINS
zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
END_2D
!
CALL iom_put( 'sig1_pnorm' , zsig1_p )
CALL iom_put( 'sig2_pnorm' , zsig2_p )
CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00 )
CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00 )
DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
......@@ -890,9 +882,6 @@ CONTAINS
! --- yieldcurve --- !
IF( iom_use('yield11') .OR. iom_use('yield12') .OR. iom_use('yield22')) THEN
CALL lbc_lnk( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp )
CALL iom_put( 'yield11', zyield11 * zmsk00 )
CALL iom_put( 'yield22', zyield22 * zmsk00 )
CALL iom_put( 'yield12', zyield12 * zmsk00 )
......@@ -900,31 +889,22 @@ CONTAINS
! --- anisotropy tensor --- !
IF( iom_use('aniso') ) THEN
CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp )
CALL iom_put( 'aniso' , paniso_11 * zmsk00 )
ENDIF
! --- SIMIP --- !
IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
& iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
!
CALL lbc_lnk( 'icedyn_rhg_eap', zspgU, 'U', -1.0_wp, zspgV, 'V', -1.0_wp, &
& zCorU, 'U', -1.0_wp, zCorV, 'V', -1.0_wp, &
& zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp )
CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x)
CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y)
CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x)
CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y)
CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x)
CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y)
ENDIF
IF( iom_use('dssh_dx') ) CALL iom_put( 'dssh_dx' , zspgU(A2D(0)) * zmsk00 ) ! Sea-surface tilt term in force balance (x)
IF( iom_use('dssh_dy') ) CALL iom_put( 'dssh_dy' , zspgV(A2D(0)) * zmsk00 ) ! Sea-surface tilt term in force balance (y)
IF( iom_use('corstrx') ) CALL iom_put( 'corstrx' , zCorU(A2D(0)) * zmsk00 ) ! Coriolis force term in force balance (x)
IF( iom_use('corstry') ) CALL iom_put( 'corstry' , zCorV(A2D(0)) * zmsk00 ) ! Coriolis force term in force balance (y)
IF( iom_use('intstrx') ) CALL iom_put( 'intstrx' , zfU (A2D(0)) * zmsk00 ) ! Internal force term in force balance (x)
IF( iom_use('intstry') ) CALL iom_put( 'intstry' , zfV (A2D(0)) * zmsk00 ) ! Internal force term in force balance (y)
IF( iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
& iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
!
ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
& zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
ALLOCATE( zdiag_xmtrp_ice(A2D(0)) , zdiag_ymtrp_ice(A2D(0)) , &
& zdiag_xmtrp_snw(A2D(0)) , zdiag_ymtrp_snw(A2D(0)) , zdiag_xatrp(A2D(0)) , zdiag_yatrp(A2D(0)) )
!
DO_2D( 0, 0, 0, 0 )
! 2D ice mass, snow mass, area transport arrays (X, Y)
......@@ -942,10 +922,6 @@ CONTAINS
END_2D
CALL lbc_lnk( 'icedyn_rhg_eap', zdiag_xmtrp_ice, 'U', -1.0_wp, zdiag_ymtrp_ice, 'V', -1.0_wp, &
& zdiag_xmtrp_snw, 'U', -1.0_wp, zdiag_ymtrp_snw, 'V', -1.0_wp, &
& zdiag_xatrp , 'U', -1.0_wp, zdiag_yatrp , 'V', -1.0_wp )
CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s)
CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport
CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s)
......@@ -962,11 +938,11 @@ CONTAINS
IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
IF( iom_use('uice_cvg') ) THEN
IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &
& ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )
CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(A2D(0)) - zu_ice(:,:) ) * zbeta(A2D(0)) * umask(A2D(0),1) , &
& ABS( v_ice(A2D(0)) - zv_ice(:,:) ) * zbeta(A2D(0)) * vmask(A2D(0),1) ) * zmsk15(:,:) )
ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
& ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(A2D(0)) - zu_ice(:,:) ) * umask(A2D(0),1) , &
& ABS( v_ice(A2D(0)) - zv_ice(:,:) ) * vmask(A2D(0),1) ) * zmsk15(:,:) )
ENDIF
ENDIF
ENDIF
......@@ -987,22 +963,27 @@ CONTAINS
!!
!! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index
REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities
REAL(wp), DIMENSION(:,:), INTENT(in) :: pmsk15
INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index
REAL(wp), DIMENSION(:,:) , INTENT(in) :: pu, pv ! now velocities
REAL(wp), DIMENSION(A2D(0)), INTENT(in) :: pub, pvb ! before velocities
REAL(wp), DIMENSION(A2D(0)), INTENT(in) :: pmsk15
!!
INTEGER :: it, idtime, istatus
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zresm ! local real
CHARACTER(len=20) :: clname
LOGICAL :: ll_maxcvg
REAL(wp), DIMENSION(A2D(0),2) :: zres
REAL(wp), DIMENSION(2) :: ztmp
!!----------------------------------------------------------------------
ll_maxcvg = .FALSE.
!
! create file
IF( kt == nit000 .AND. kiter == 1 ) THEN
!
IF( lwp ) THEN
WRITE(numout,*)
WRITE(numout,*) 'rhg_cvg_eap : ice rheology convergence control'
WRITE(numout,*) 'rhg_cvg : ice rheology convergence control'
WRITE(numout,*) '~~~~~~~'
ENDIF
!
......@@ -1011,7 +992,7 @@ CONTAINS
IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime )
istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid )
istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid )
istatus = NF90_ENDDEF(ncvgid)
ENDIF
!
......@@ -1025,11 +1006,21 @@ CONTAINS
zresm = 0._wp
ELSE
zresm = 0._wp
DO_2D( 0, 0, 0, 0 )
zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
& ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) )
END_2D
CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain
IF( ll_maxcvg ) THEN ! error max over the domain
DO_2D( 0, 0, 0, 0 )
zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
& ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) )
END_2D
CALL mpp_max( 'icedyn_rhg_eap', zresm )
ELSE ! error averaged over the domain
DO_2D( 0, 0, 0, 0 )
zres(ji,jj,1) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
& ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj)
zres(ji,jj,2) = pmsk15(ji,jj)
END_2D
ztmp(:) = glob_sum_vec( 'icedyn_rhg_eap', zres )
IF( ztmp(2) /= 0._wp ) zresm = ztmp(1) / ztmp(2)
ENDIF
ENDIF
IF( lwm ) THEN
......
src/ICE/icedyn_rhg_evp.F90
View file @ 3efbcb17
......@@ -114,10 +114,10 @@ CONTAINS
!! Bouillon et al., Ocean Modelling 2013
!! Kimmritz et al., Ocean Modelling 2016 & 2017
!!-------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i !
REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i !
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i !
REAL(wp), DIMENSION(A2D(0)), INTENT( out) :: pshear_i , pdivu_i , pdelta_i !
!!
INTEGER :: ji, jj ! dummy loop indices
INTEGER :: jter ! local integers
......@@ -135,38 +135,39 @@ CONTAINS
!
REAL(wp) :: zfac_x, zfac_y
!
REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points
REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017
REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta, P/delta at T points
REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017
!
REAL(wp), DIMENSION(jpi,jpj) :: zdt_m ! (dt / ice-snow_mass) on T points
REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points
REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points
REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points
REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points
REAL(wp), DIMENSION(jpi,jpj) :: zdt_m ! (dt / ice-snow_mass) on T points
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: zaU , zaV ! ice fraction on U/V points
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points
REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points
!
REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear
REAL(wp), DIMENSION(jpi,jpj) :: zten_i, zshear ! tension, shear
REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components
REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope:
! ! ocean surface (ssh_m) if ice is not embedded
! ! ice bottom surface if ice is embedded
REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses
REAL(wp), DIMENSION(jpi,jpj) :: zspgU, zspgV ! surface pressure gradient at U/V points
REAL(wp), DIMENSION(jpi,jpj) :: zCorU, zCorV ! Coriolis stress array
REAL(wp), DIMENSION(jpi,jpj) :: ztaux_ai, ztauy_ai ! ice-atm. stress at U-V points
REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! ice-ocean stress at U-V points
REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast)
REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast)
REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear
REAL(wp), DIMENSION(A2D(0)) :: zten_i, zshear ! tension, shear
REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components
REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope:
! ! ocean surface (ssh_m) if ice is not embedded
! ! ice bottom surface if ice is embedded
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: zfU , zfV ! internal stresses
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: zspgU, zspgV ! surface pressure gradient at U/V points
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: zCorU, zCorV ! Coriolis stress array
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: ztaux_ai, ztauy_ai ! ice-atm. stress at U-V points
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: ztaux_oi, ztauy_oi ! ice-ocean stress at U-V points
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast)
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast)
!
REAL(wp), DIMENSION(jpi,jpj) :: zmsk, zmsk00, zmsk15
REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays
REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence
REAL(wp), DIMENSION(jpi,jpj) :: zmsk
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: zmsk01x, zmsk01y ! dummy arrays
REAL(wp), DIMENSION(A2D(nn_hls-1)) :: zmsk00x, zmsk00y ! mask for ice presence
REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter
REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small
REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small
!! --- check convergence
REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice
REAL(wp), DIMENSION(A2D(0)) :: zmsk00, zmsk15
REAL(wp), DIMENSION(A2D(0)) :: zu_ice, zv_ice
!! --- diags
REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p
......@@ -186,13 +187,15 @@ CONTAINS
IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
!
! for diagnostics and convergence tests
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice
zmsk (ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 1 if ice , 0 if no ice
zmsk (ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 1 if ice , 0 if no ice
END_2D
! for diagnostics and convergence tests
DO_2D( 0, 0, 0, 0 )
zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice
END_2D
IF( nn_rhg_chkcvg > 0 ) THEN
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
END_2D
ENDIF
......@@ -265,6 +268,7 @@ CONTAINS
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zm1 = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) ) ! Ice/snow mass at U-V points
!!$ zm1 = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) + rhow * (vt_ip(ji,jj) + vt_il(ji,jj)) ) ! clem: this should replace the above
zmf (ji,jj) = zm1 * ff_t(ji,jj) ! Coriolis at T points (m*f)
zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) ! dt/m at T points (for alpha and beta coefficients)
END_2D
......@@ -277,8 +281,11 @@ CONTAINS
! Ice/snow mass at U-V points
zm1 = ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ) )
!!$ zm1 = ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ) + rhow * (vt_ip(ji ,jj ) + vt_il(ji ,jj )) ) ! clem: this should replace the above
zm2 = ( rhos * vt_s(ji+1,jj ) + rhoi * vt_i(ji+1,jj ) )
!!$ zm2 = ( rhos * vt_s(ji+1,jj ) + rhoi * vt_i(ji+1,jj ) + rhow * (vt_ip(ji+1,jj ) + vt_il(ji+1,jj )) ) ! clem: this should replace the above
zm3 = ( rhos * vt_s(ji ,jj+1) + rhoi * vt_i(ji ,jj+1) )
!!$ zm3 = ( rhos * vt_s(ji ,jj+1) + rhoi * vt_i(ji ,jj+1) + rhow * (vt_ip(ji ,jj+1) + vt_il(ji ,jj+1)) ) ! clem: this should replace the above
zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
......@@ -292,8 +299,12 @@ CONTAINS
zmV_t(ji,jj) = zmassV * z1_dtevp
! Drag ice-atm.
ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
! and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
ztaux_ai(ji,jj) = zaU(ji,jj) * 0.5_wp * ( utau_ice(ji,jj) + utau_ice(ji+1,jj) ) * &
& ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
ztauy_ai(ji,jj) = zaV(ji,jj) * 0.5_wp * ( vtau_ice(ji,jj) + vtau_ice(ji,jj+1) ) * &
& ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
......@@ -324,11 +335,12 @@ CONTAINS
! ice-bottom stress at V points
zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
END_2D
DO_2D( 0, 0, 0, 0 )
! ice_bottom stress at T points
zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
END_2D
CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1.0_wp )
!
ELSE !-- no landfast
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
......@@ -348,7 +360,7 @@ CONTAINS
!
! convergence test
IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
END_2D
......@@ -718,7 +730,7 @@ CONTAINS
! ! ==================== !
END DO ! end loop over jter !
! ! ==================== !
IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta )
IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta(A2D(0)) )
!
IF( ll_advups .AND. ln_str_H79 ) CALL lbc_lnk( 'icedyn_rhg_evp', strength, 'T', 1.0_wp )
!
......@@ -772,8 +784,7 @@ CONTAINS
END_2D
CALL lbc_lnk( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, &
& zshear , 'T', 1._wp, zdelta , 'T', 1._wp, zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12, 'F', 1._wp )
CALL lbc_lnk( 'icedyn_rhg_evp', zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp )
! --- Store the stress tensor for the next time step --- !
pstress1_i (:,:) = zs1 (:,:)
......@@ -783,34 +794,25 @@ CONTAINS
! 5) diagnostics
!------------------------------------------------------------------------------!
! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
& iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
!
CALL lbc_lnk( 'icedyn_rhg_evp', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, &
& ztaux_ai, 'U', -1.0_wp, ztauy_ai, 'V', -1.0_wp, &
& ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp )
!
CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
ENDIF
IF( iom_use('utau_oi') ) CALL iom_put( 'utau_oi' , ztaux_oi(A2D(0)) * zmsk00 )
IF( iom_use('vtau_oi') ) CALL iom_put( 'vtau_oi' , ztauy_oi(A2D(0)) * zmsk00 )
IF( iom_use('utau_ai') ) CALL iom_put( 'utau_ai' , ztaux_ai(A2D(0)) * zmsk00 )
IF( iom_use('vtau_ai') ) CALL iom_put( 'vtau_ai' , ztauy_ai(A2D(0)) * zmsk00 )
IF( iom_use('utau_bi') ) CALL iom_put( 'utau_bi' , ztaux_bi(A2D(0)) * zmsk00 )
IF( iom_use('vtau_bi') ) CALL iom_put( 'vtau_bi' , ztauy_bi(A2D(0)) * zmsk00 )
! --- divergence, shear and strength --- !
IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence
IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear
IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength
IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , zdelta * zmsk00 ) ! delta
IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i (A2D(0)) * zmsk00 ) ! divergence
IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i(A2D(0)) * zmsk00 ) ! shear
IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength(A2D(0)) * zmsk00 ) ! strength
IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , zdelta (A2D(0)) * zmsk00 ) ! delta
! --- Stress tensor invariants (SIMIP diags) --- !
IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
!
ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
ALLOCATE( zsig_I(A2D(0)) , zsig_II(A2D(0)) )
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
! Ice stresses
! sigma1, sigma2, sigma12 are some recombination of the stresses (HD MWR002, Bouillon et al., OM2013)
! not to be confused with stress tensor components, stress invariants, or stress principal components
......@@ -825,8 +827,8 @@ CONTAINS
END_2D
!
IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress
IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00 ) ! Normal stress
IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00 ) ! Maximum shear stress
DEALLOCATE ( zsig_I, zsig_II )
......@@ -838,9 +840,9 @@ CONTAINS
! Recommendation 2 : for EVP, no need to use viscosities at last iteration (stress is properly iterated)
IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
!
ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
ALLOCATE( zsig1_p(A2D(0)) , zsig2_p(A2D(0)) , zsig_I(A2D(0)) , zsig_II(A2D(0)) )
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
! For EVP solvers, ice stresses at current iterates can be used
! following Lemieux & Dupont (2020)
......@@ -859,33 +861,26 @@ CONTAINS
zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
END_2D
!
CALL iom_put( 'sig1_pnorm' , zsig1_p * zmsk00 )
CALL iom_put( 'sig2_pnorm' , zsig2_p * zmsk00 )
CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00 )
CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00 )
DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
ENDIF
! --- SIMIP --- !
IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
& iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
!
CALL lbc_lnk( 'icedyn_rhg_evp', zspgU, 'U', -1.0_wp, zspgV, 'V', -1.0_wp, &
& zCorU, 'U', -1.0_wp, zCorV, 'V', -1.0_wp, zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp )
CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x)
CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y)
CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x)
CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y)
CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x)
CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y)
ENDIF
IF( iom_use('dssh_dx') ) CALL iom_put( 'dssh_dx' , zspgU(A2D(0)) * zmsk00 ) ! Sea-surface tilt term in force balance (x)
IF( iom_use('dssh_dy') ) CALL iom_put( 'dssh_dy' , zspgV(A2D(0)) * zmsk00 ) ! Sea-surface tilt term in force balance (y)
IF( iom_use('corstrx') ) CALL iom_put( 'corstrx' , zCorU(A2D(0)) * zmsk00 ) ! Coriolis force term in force balance (x)
IF( iom_use('corstry') ) CALL iom_put( 'corstry' , zCorV(A2D(0)) * zmsk00 ) ! Coriolis force term in force balance (y)
IF( iom_use('intstrx') ) CALL iom_put( 'intstrx' , zfU (A2D(0)) * zmsk00 ) ! Internal force term in force balance (x)
IF( iom_use('intstry') ) CALL iom_put( 'intstry' , zfV (A2D(0)) * zmsk00 ) ! Internal force term in force balance (y)
IF( iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
& iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
!
ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
& zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
ALLOCATE( zdiag_xmtrp_ice(A2D(0)) , zdiag_ymtrp_ice(A2D(0)) , &
& zdiag_xmtrp_snw(A2D(0)) , zdiag_ymtrp_snw(A2D(0)) , zdiag_xatrp(A2D(0)) , zdiag_yatrp(A2D(0)) )
!
DO_2D( 0, 0, 0, 0 )
! 2D ice mass, snow mass, area transport arrays (X, Y)
......@@ -903,10 +898,6 @@ CONTAINS
END_2D
CALL lbc_lnk( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1.0_wp, zdiag_ymtrp_ice, 'V', -1.0_wp, &
& zdiag_xmtrp_snw, 'U', -1.0_wp, zdiag_ymtrp_snw, 'V', -1.0_wp, &
& zdiag_xatrp , 'U', -1.0_wp, zdiag_yatrp , 'V', -1.0_wp )
CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s)
CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport
CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s)
......@@ -923,11 +914,11 @@ CONTAINS
IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
IF( iom_use('uice_cvg') ) THEN
IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &
& ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )
CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(A2D(0)) - zu_ice(:,:) ) * zbeta(A2D(0)) * umask(A2D(0),1) , &
& ABS( v_ice(A2D(0)) - zv_ice(:,:) ) * zbeta(A2D(0)) * vmask(A2D(0),1) ) * zmsk15(:,:) )
ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
& ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(A2D(0)) - zu_ice(:,:) ) * umask(A2D(0),1) , &
& ABS( v_ice(A2D(0)) - zv_ice(:,:) ) * vmask(A2D(0),1) ) * zmsk15(:,:) )
ENDIF
ENDIF
ENDIF
......@@ -948,17 +939,18 @@ CONTAINS
!!
!! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index
REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities
REAL(wp), DIMENSION(:,:), INTENT(in) :: pmsk15
INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index
REAL(wp), DIMENSION(:,:) , INTENT(in) :: pu, pv ! now velocities
REAL(wp), DIMENSION(A2D(0)), INTENT(in) :: pub, pvb ! before velocities
REAL(wp), DIMENSION(A2D(0)), INTENT(in) :: pmsk15
!!
INTEGER :: it, idtime, istatus
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zresm ! local real
CHARACTER(len=20) :: clname
LOGICAL :: ll_maxcvg
REAL(wp), DIMENSION(jpi,jpj,2) :: zres
REAL(wp), DIMENSION(2) :: ztmp
REAL(wp), DIMENSION(A2D(0),2) :: zres
REAL(wp), DIMENSION(2) :: ztmp
!!----------------------------------------------------------------------
ll_maxcvg = .FALSE.
!
......
src/ICE/icedyn_rhg_vp.F90
View file @ 3efbcb17
......@@ -332,7 +332,10 @@ CONTAINS
v_oceU(ji,jj) = 0.25_wp * ( (v_oce(ji,jj) + v_oce(ji,jj-1)) + (v_oce(ji+1,jj) + v_oce(ji+1,jj-1)) ) * umask(ji,jj,1)
! Wind stress
ztaux_ai(ji,jj) = za_iU(ji,jj) * utau_ice(ji,jj)
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
! and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
ztaux_ai(ji,jj) = za_iU(ji,jj) * 0.5_wp * ( utau_ice(ji,jj) + utau_ice(ji+1,jj) ) * &
& ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
! Force due to sea surface tilt(- m*g*GRAD(ssh))
zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
......@@ -369,7 +372,10 @@ CONTAINS
u_oceV(ji,jj) = 0.25_wp * ( (u_oce(ji,jj) + u_oce(ji-1,jj)) + (u_oce(ji,jj+1) + u_oce(ji-1,jj+1)) ) * vmask(ji,jj,1)
! Wind stress
ztauy_ai(ji,jj) = za_iV(ji,jj) * vtau_ice(ji,jj)
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
! and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
ztauy_ai(ji,jj) = za_iV(ji,jj) * 0.5_wp * ( vtau_ice(ji,jj) + vtau_ice(ji,jj+1) ) * &
& ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
! Force due to sea surface tilt(- m*g*GRAD(ssh))
zspgV(ji,jj) = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
......
src/ICE/iceistate.F90
View file @ 3efbcb17
......@@ -37,6 +37,7 @@ MODULE iceistate
USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE fldread ! read input fields
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
# if defined key_agrif
USE agrif_oce
......@@ -113,11 +114,11 @@ CONTAINS
INTEGER :: ji, jj, jk, jl ! dummy loop indices
REAL(wp) :: ztmelts, zsshadj, area
INTEGER , DIMENSION(4) :: itest
REAL(wp), DIMENSION(jpi,jpj) :: z2d
REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator
REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, ztm_s_ini !data from namelist or nc file
REAL(wp), DIMENSION(jpi,jpj) :: zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
REAL(wp), DIMENSION(jpi,jpj) :: zapnd_ini, zhpnd_ini, zhlid_ini !data from namelist or nc file
REAL(wp), DIMENSION(A2D(0)) :: zmsk ! ice indicator
REAL(wp), DIMENSION(A2D(0)) :: zht_i_ini, zat_i_ini, ztm_s_ini !data from namelist or nc file
REAL(wp), DIMENSION(A2D(0)) :: zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
REAL(wp), DIMENSION(A2D(0)) :: zapnd_ini, zhpnd_ini, zhlid_ini !data from namelist or nc file
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zti_3d , zts_3d !temporary arrays
!!
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d
......@@ -135,53 +136,59 @@ CONTAINS
CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)
!
! surface temperature and conductivity
DO jl = 1, jpl
t_su (:,:,jl) = rt0 * tmask(:,:,1) ! temp at the surface
cnd_ice(:,:,jl) = 0._wp ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
END DO
!
! ice and snw temperatures
DO jl = 1, jpl
DO jk = 1, nlay_i
t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
END DO
DO jk = 1, nlay_s
t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
END DO
END DO
!
! specific temperatures for coupled runs
tn_ice (:,:,:) = t_i (:,:,1,:)
t1_ice (:,:,:) = t_i (:,:,1,:)
! heat contents
e_i (:,:,:,:) = 0._wp
e_s (:,:,:,:) = 0._wp
! general fields
a_i (:,:,:) = 0._wp
v_i (:,:,:) = 0._wp
v_s (:,:,:) = 0._wp
sv_i(:,:,:) = 0._wp
oa_i(:,:,:) = 0._wp
!
h_i (:,:,:) = 0._wp
h_s (:,:,:) = 0._wp
s_i (:,:,:) = 0._wp
o_i (:,:,:) = 0._wp
!
! melt ponds
a_ip (:,:,:) = 0._wp
v_ip (:,:,:) = 0._wp
v_il (:,:,:) = 0._wp
a_ip_eff (:,:,:) = 0._wp
h_ip (:,:,:) = 0._wp
h_il (:,:,:) = 0._wp
! == reduced arrays == !
DO_2D( 0, 0, 0, 0 )
!
cnd_ice(ji,jj,jl) = 0._wp ! conductivity at the ice top
!
tn_ice(ji,jj,jl) = t_i (ji,jj,1,jl) ! temp for coupled runs
t1_ice(ji,jj,jl) = t_i (ji,jj,1,jl) ! temp for coupled runs
!
a_ip_eff(ji,jj,jl) = 0._wp ! melt pond effective fraction
END_2D
!
! == full arrays == !
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
! heat contents
DO jk = 1, nlay_i
e_i(ji,jj,jk,jl) = 0._wp
t_i(ji,jj,jk,jl) = rt0 * tmask(ji,jj,1) ! ice temp
ENDDO
DO jk = 1, nlay_s
e_s(ji,jj,jk,jl) = 0._wp
t_s(ji,jj,jk,jl) = rt0 * tmask(ji,jj,1) ! snw temp
ENDDO
!
! general fields
a_i (ji,jj,jl) = 0._wp
v_i (ji,jj,jl) = 0._wp
v_s (ji,jj,jl) = 0._wp
sv_i(ji,jj,jl) = 0._wp
oa_i(ji,jj,jl) = 0._wp
h_i (ji,jj,jl) = 0._wp
h_s (ji,jj,jl) = 0._wp
s_i (ji,jj,jl) = 0._wp
o_i (ji,jj,jl) = 0._wp
t_su(ji,jj,jl) = rt0 * tmask(ji,jj,1)
!
! melt ponds
a_ip(ji,jj,jl) = 0._wp
v_ip(ji,jj,jl) = 0._wp
v_il(ji,jj,jl) = 0._wp
h_ip(ji,jj,jl) = 0._wp
h_il(ji,jj,jl) = 0._wp
!
END_2D
!
ENDDO
!
! ice velocities
u_ice (:,:) = 0._wp
v_ice (:,:) = 0._wp
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
u_ice(ji,jj) = 0._wp
v_ice(ji,jj) = 0._wp
END_2D
!
!------------------------------------------------------------------------
! 2) overwrite some of the fields with namelist parameters or netcdf file
......@@ -194,30 +201,30 @@ CONTAINS
! !---------------!
IF( nn_iceini_file == 1 )THEN ! Read a file !
! !---------------!
WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp
ELSEWHERE ; zswitch(:,:) = 0._wp
WHERE( ff_t(A2D(0)) >= 0._wp ) ; zmsk(:,:) = 1._wp
ELSEWHERE ; zmsk(:,:) = 0._wp
END WHERE
!
CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
!
! -- mandatory fields -- !
zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1)
zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1)
zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1)
zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * smask0(:,:)
zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * smask0(:,:)
zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * smask0(:,:)
! -- optional fields -- !
! if fields do not exist then set them to the values present in the namelist (except for temperatures)
!
! ice salinity
IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
& si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
& si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zmsk + rn_smi_ini_s * (1._wp - zmsk) ) * smask0(:,:)
!
! temperatures
IF ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. &
& TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN
si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zmsk + rn_tmi_ini_s * (1._wp - zmsk) ) * smask0(:,:)
si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zmsk + rn_tsu_ini_s * (1._wp - zmsk) ) * smask0(:,:)
si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zmsk + rn_tms_ini_s * (1._wp - zmsk) ) * smask0(:,:)
ENDIF
IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
& si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 )
......@@ -234,67 +241,60 @@ CONTAINS
!
! pond concentration
IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
& si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc.
& si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zmsk + rn_apd_ini_s * (1._wp - zmsk) ) * smask0(:,:) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc.
& * si(jp_ati)%fnow(:,:,1)
!
! pond depth
IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
& si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
& si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zmsk + rn_hpd_ini_s * (1._wp - zmsk) ) * smask0(:,:)
!
! pond lid depth
IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) &
& si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zswitch + rn_hld_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
& si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zmsk + rn_hld_ini_s * (1._wp - zmsk) ) * smask0(:,:)
!
zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1)
ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1)
zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1)
ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1)
zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1)
zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1)
zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1) * tmask(:,:,1)
zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * smask0(:,:)
ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * smask0(:,:)
zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * smask0(:,:)
ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * smask0(:,:)
zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * smask0(:,:)
zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * smask0(:,:)
zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1) * smask0(:,:)
!
! change the switch for the following
WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
ELSEWHERE ; zswitch(:,:) = 0._wp
END WHERE
! !---------------!
ELSE ! Read namelist !
! !---------------!
! no ice if (sst - Tfreez) >= thresold
WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp
ELSEWHERE ; zswitch(:,:) = tmask(:,:,1)
WHERE( ( sst_m(A2D(0)) - (t_bo(A2D(0)) - rt0) ) * smask0(:,:) >= rn_thres_sst ) ; zmsk(:,:) = 0._wp
ELSEWHERE ; zmsk(:,:) = smask0(:,:)
END WHERE
!
! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
WHERE( ff_t(:,:) >= 0._wp )
zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:)
WHERE( ff_t(A2D(0)) >= 0._wp )
zht_i_ini(:,:) = rn_hti_ini_n * zmsk(:,:)
zht_s_ini(:,:) = rn_hts_ini_n * zmsk(:,:)
zat_i_ini(:,:) = rn_ati_ini_n * zmsk(:,:)
zsm_i_ini(:,:) = rn_smi_ini_n * zmsk(:,:)
ztm_i_ini(:,:) = rn_tmi_ini_n * zmsk(:,:)
zt_su_ini(:,:) = rn_tsu_ini_n * zmsk(:,:)
ztm_s_ini(:,:) = rn_tms_ini_n * zmsk(:,:)
zapnd_ini(:,:) = rn_apd_ini_n * zmsk(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
zhpnd_ini(:,:) = rn_hpd_ini_n * zmsk(:,:)
zhlid_ini(:,:) = rn_hld_ini_n * zmsk(:,:)
ELSEWHERE
zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:)
zht_i_ini(:,:) = rn_hti_ini_s * zmsk(:,:)
zht_s_ini(:,:) = rn_hts_ini_s * zmsk(:,:)
zat_i_ini(:,:) = rn_ati_ini_s * zmsk(:,:)
zsm_i_ini(:,:) = rn_smi_ini_s * zmsk(:,:)
ztm_i_ini(:,:) = rn_tmi_ini_s * zmsk(:,:)
zt_su_ini(:,:) = rn_tsu_ini_s * zmsk(:,:)
ztm_s_ini(:,:) = rn_tms_ini_s * zmsk(:,:)
zapnd_ini(:,:) = rn_apd_ini_s * zmsk(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
zhpnd_ini(:,:) = rn_hpd_ini_s * zmsk(:,:)
zhlid_ini(:,:) = rn_hld_ini_s * zmsk(:,:)
END WHERE
!
ENDIF
! make sure ponds = 0 if no ponds scheme
IF ( .NOT.ln_pnd ) THEN
zapnd_ini(:,:) = 0._wp
......@@ -311,7 +311,7 @@ CONTAINS
!----------------!
! select ice covered grid points
npti = 0 ; nptidx(:) = 0
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
npti = npti + 1
nptidx(npti) = (jj - 1) * jpi + ji
......@@ -363,6 +363,16 @@ CONTAINS
DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
& zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d )
! this call is needed because of the calculations above that are done only in the interior
CALL lbc_lnk( 'iceistate', a_i , 'T', 1._wp, h_i , 'T', 1._wp, h_s , 'T', 1._wp, &
& zti_3d, 'T', 1._wp, zts_3d, 'T', 1._wp, t_su, 'T', 1._wp, &
& s_i , 'T', 1._wp, a_ip , 'T', 1._wp, h_ip, 'T', 1._wp, h_il, 'T', 1._wp )
! switch for the following
WHERE( SUM(a_i(:,:,:),dim=3) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
ELSEWHERE ; zswitch(:,:) = 0._wp
END WHERE
! calculate extensive and intensive variables
CALL ice_var_salprof ! for sz_i
DO jl = 1, jpl
......@@ -398,15 +408,17 @@ CONTAINS
ENDIF
#endif
! Melt ponds
WHERE( a_i > epsi10 ) ; a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
ELSEWHERE ; a_ip_eff(:,:,:) = 0._wp
WHERE( a_i(A2D(0),:) > epsi10 ) ; a_ip_eff(:,:,:) = a_ip(A2D(0),:) / a_i(A2D(0),:)
ELSEWHERE ; a_ip_eff(:,:,:) = 0._wp
END WHERE
v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:)
! specific temperatures for coupled runs
tn_ice(:,:,:) = t_su(:,:,:)
t1_ice(:,:,:) = t_i (:,:,1,:)
DO_2D( 0, 0, 0, 0 )
tn_ice(ji,jj,:) = t_su(ji,jj,:)
t1_ice(ji,jj,:) = t_i (ji,jj,1,:)
END_2D
!
! ice concentration should not exceed amax
at_i(:,:) = SUM( a_i, dim=3 )
......@@ -546,8 +558,8 @@ CONTAINS
ENDIF
!
DO ifpr = 1, jpfldi
ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
ALLOCATE( si(ifpr)%fnow(A2D(0),1) )
IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(A2D(0),1,2) )
END DO
!
! fill si with slf_i and control print
......
src/ICE/iceitd.F90
View file @ 3efbcb17
......@@ -29,6 +29,7 @@ MODULE iceitd
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE prtctl ! Print control
USE timing ! Timing
USE lbclnk ! lateral boundary conditions (or mpp links)
IMPLICIT NONE
PRIVATE
......@@ -100,7 +101,7 @@ CONTAINS
at_i(:,:) = SUM( a_i, dim=3 )
!
npti = 0 ; nptidx(:) = 0
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF ( at_i(ji,jj) > epsi10 ) THEN
npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji
......@@ -327,6 +328,9 @@ CONTAINS
!
ENDIF
!
! the following fields need to be updated in the halos (done afterwards):
! a_i, v_i, v_s, sv_i, oa_i, h_i, a_ip, v_ip, v_il, t_su, e_i, e_s
!
IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_rem', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft)
IF( ln_timing ) CALL timing_stop ('iceitd_rem')
......@@ -626,7 +630,7 @@ CONTAINS
DO jl = 1, jpl-1 ! identify thicknesses that are too big
! !---------------------------------------
npti = 0 ; nptidx(:) = 0
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF( a_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji
......@@ -662,7 +666,7 @@ CONTAINS
DO jl = jpl-1, 1, -1 ! Identify thicknesses that are too small
! !-----------------------------------------
npti = 0 ; nptidx(:) = 0
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF( a_i(ji,jj,jl+1) > 0._wp .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN
npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji
......@@ -687,6 +691,14 @@ CONTAINS
!
END DO
!
! clem: the 2 lbc below could be avoided if calculations above were performed over the full domain
! => decide if it is more efficient this way or not?
! note: ice_itd_reb is called in icedyn
! and in icethd (but once the arrays are already updated on the boundaries)
CALL lbc_lnk( 'iceitd_reb', a_i, 'T', 1._wp, v_i , 'T', 1._wp, v_s , 'T', 1._wp, sv_i, 'T', 1._wp, oa_i, 'T', 1._wp, &
& h_i, 'T', 1._wp, a_ip, 'T', 1._wp, v_ip, 'T', 1._wp, v_il, 'T', 1._wp, t_su, 'T', 1._wp )
CALL lbc_lnk( 'iceitd_reb', e_i, 'T', 1._wp, e_s , 'T', 1._wp )
!
IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_reb', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft)
IF( ln_timing ) CALL timing_stop ('iceitd_reb')
......
src/ICE/icesbc.F90
View file @ 3efbcb17
......@@ -58,7 +58,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: utau_ice, vtau_ice ! air-ice stress [N/m2]
!!
INTEGER :: ji, jj ! dummy loop index
REAL(wp), DIMENSION(jpi,jpj) :: zutau_ice, zvtau_ice
REAL(wp), DIMENSION(A2D(0)) :: zutau_ice, zvtau_ice
!!-------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('icesbc')
......@@ -70,25 +70,38 @@ CONTAINS
ENDIF
!
SELECT CASE( ksbc )
CASE( jp_usr ) ; CALL usrdef_sbc_ice_tau( kt ) ! user defined formulation
CASE( jp_blk )
CALL blk_ice_1( sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), &
& theta_air_zt(:,:), q_air_zt(:,:), & ! #LB: known from "sbc_oce" module...
& sf(jp_slp )%fnow(:,:,1), u_ice, v_ice, tm_su , & ! inputs
& putaui = utau_ice, pvtaui = vtau_ice ) ! outputs
! CASE( jp_abl ) utau_ice & vtau_ice are computed in ablmod
CASE( jp_purecpl ) ; CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) ! Coupled formulation
!
CASE( jp_usr ) !--- User defined formulation
!
CALL usrdef_sbc_ice_tau( kt )
!
CASE( jp_blk ) !--- Forced formulation
!
CALL blk_ice_1( sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), theta_air_zt(:,:), q_air_zt(:,:), & ! <<== in
& sf(jp_slp )%fnow(:,:,1), u_ice(A2D(0)), v_ice(A2D(0)), tm_su(:,:), & ! <<== in
& putaui = utau_ice(A2D(0)), pvtaui = vtau_ice(A2D(0)) ) ! ==>> out
!
!CASE( jp_abl ) !--- ABL formulation (utau_ice & vtau_ice are computed in ablmod)
!
CASE( jp_purecpl ) !--- Coupled formulation
!
CALL sbc_cpl_ice_tau( utau_ice(A2D(0)) , vtau_ice(A2D(0)) )
!
END SELECT
!
IF( ln_mixcpl) THEN ! Case of a mixed Bulk/Coupled formulation
CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
IF( ln_mixcpl) THEN !--- Case of a mixed Bulk/Coupled formulation
!
CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
!
DO_2D( 0, 0, 0, 0 )
utau_ice(ji,jj) = utau_ice(ji,jj) * xcplmask(ji,jj,0) + zutau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
END_2D
CALL lbc_lnk( 'icesbc', utau_ice, 'U', -1.0_wp, vtau_ice, 'V', -1.0_wp )
!
ENDIF
!
CALL lbc_lnk( 'icesbc', utau_ice, 'U', -1.0_wp, vtau_ice, 'V', -1.0_wp )
!
IF( ln_timing ) CALL timing_stop('icesbc')
!
END SUBROUTINE ice_sbc_tau
......@@ -129,25 +142,49 @@ CONTAINS
WRITE(numout,*)'~~~~~~~~~~~~~~~'
ENDIF
! !== ice albedo ==!
CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, cloud_fra, alb_ice )
CALL ice_alb( ln_pnd_alb, t_su(A2D(0),:), h_i(A2D(0),:), h_s(A2D(0),:), a_ip_eff(:,:,:), h_ip(A2D(0),:), cloud_fra(:,:), & ! <<== in
& alb_ice(:,:,:) ) ! ==>> out
!
SELECT CASE( ksbc ) !== fluxes over sea ice ==!
!
CASE( jp_usr ) !--- user defined formulation
!
CALL usrdef_sbc_ice_flx( kt, h_s, h_i )
!
CASE( jp_blk, jp_abl ) !--- bulk formulation & ABL formulation
CALL blk_ice_2 ( t_su, h_s, h_i, alb_ice, &
& theta_air_zt(:,:), q_air_zt(:,:), & ! #LB: known from "sbc_oce" module...
& sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), &
& sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) )
IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( kt, picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
IF( nn_flxdist /= -1 ) CALL ice_flx_dist ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
!
CALL blk_ice_2( t_su(A2D(0),:), h_s(A2D(0),:), h_i(A2D(0),:), & ! <<== in
& alb_ice(:,:,:), theta_air_zt(:,:), q_air_zt(:,:), & ! <<== in
& sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), & ! <<== in
& sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) ) ! <<== in
!
IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( kt, picefr=at_i_b(:,:), palbi=alb_ice(:,:,:), &
& psst=sst_m(A2D(0)), pist=t_su(A2D(0),:), &
& phs=h_s(A2D(0),:), phi=h_i(A2D(0),:) )
CALL lbc_lnk( 'icesbc', t_su, 'T', 1.0_wp ) ! clem: t_su is needed for Met-Office only => necessary?
!
IF( nn_flxdist /= -1 ) CALL ice_flx_dist( nn_flxdist, at_i(A2D(0)), a_i(A2D(0),:), t_su(A2D(0),:), alb_ice(:,:,:), & ! <<== in
& qns_ice(:,:,:), qsr_ice(:,:,:), dqns_ice(:,:,:), & ! ==>> inout
& evap_ice(:,:,:), devap_ice(:,:,:) ) ! ==>> inout
!
! ! compute conduction flux and surface temperature (as in Jules surface module)
IF( ln_cndflx .AND. .NOT.ln_cndemulate ) &
& CALL blk_ice_qcn ( ln_virtual_itd, t_su, t_bo, h_s, h_i )
IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN
CALL blk_ice_qcn( ln_virtual_itd, t_bo(A2D(0)), h_s(A2D(0),:), h_i(A2D(0),:), & ! <<== in
& qcn_ice(:,:,:), qml_ice(:,:,:), & ! ==>> out
& qns_ice(:,:,:), t_su(A2D(0),:) ) ! ==>> inout
CALL lbc_lnk( 'icesbc', t_su, 'T', 1.0_wp ) ! clem: t_su is updated in ice_qcn => necessary?
ENDIF
!
CASE ( jp_purecpl ) !--- coupled formulation
CALL sbc_cpl_ice_flx( kt, picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
IF( nn_flxdist /= -1 ) CALL ice_flx_dist ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
!
CALL sbc_cpl_ice_flx( kt, picefr=at_i_b(:,:), palbi=alb_ice(:,:,:), psst=sst_m(A2D(0)), &
& pist=t_su(A2D(0),:), phs=h_s(A2D(0),:), phi=h_i(A2D(0),:) )
CALL lbc_lnk( 'icesbc', t_su, 'T', 1.0_wp ) ! clem: t_su is needed for Met-Office only => necessary?
!
IF( nn_flxdist /= -1 ) CALL ice_flx_dist( nn_flxdist, at_i(A2D(0)), a_i(A2D(0),:), t_su(A2D(0),:), alb_ice(:,:,:), & ! <<== in
& qns_ice(:,:,:), qsr_ice(:,:,:), dqns_ice(:,:,:), & ! ==>> inout
& evap_ice(:,:,:), devap_ice(:,:,:) ) ! ==>> inout
!
END SELECT
! !== some fluxes at the ice-ocean interface and in the leads
CALL ice_flx_other
......@@ -157,7 +194,8 @@ CONTAINS
END SUBROUTINE ice_sbc_flx
SUBROUTINE ice_flx_dist( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_flxdist )
SUBROUTINE ice_flx_dist( k_flxdist, pat_i, pa_i, ptn_ice, palb_ice, &
& pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_flx_dist ***
!!
......@@ -173,18 +211,20 @@ CONTAINS
!! using T-ice and albedo sensitivity
!! = 2 Redistribute a single flux over categories
!!-------------------------------------------------------------------
INTEGER , INTENT(in ) :: k_flxdist ! redistributor
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptn_ice ! ice surface temperature
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb_ice ! ice albedo
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqns_ice ! non solar flux
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqsr_ice ! net solar flux
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdqn_ice ! non solar flux sensitivity
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pevap_ice ! sublimation
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdevap_ice ! sublimation sensitivity
INTEGER , INTENT(in ) :: k_flxdist ! redistributor
REAL(wp), DIMENSION(A2D(0)) , INTENT(in ) :: pat_i ! ice concentration
REAL(wp), DIMENSION(A2D(0),jpl), INTENT(in ) :: pa_i ! ice concentration
REAL(wp), DIMENSION(A2D(0),jpl), INTENT(in ) :: ptn_ice ! ice surface temperature
REAL(wp), DIMENSION(A2D(0),jpl), INTENT(in ) :: palb_ice ! ice albedo
REAL(wp), DIMENSION(A2D(0),jpl), INTENT(inout) :: pqns_ice ! non solar flux
REAL(wp), DIMENSION(A2D(0),jpl), INTENT(inout) :: pqsr_ice ! net solar flux
REAL(wp), DIMENSION(A2D(0),jpl), INTENT(inout) :: pdqn_ice ! non solar flux sensitivity
REAL(wp), DIMENSION(A2D(0),jpl), INTENT(inout) :: pevap_ice ! sublimation
REAL(wp), DIMENSION(A2D(0),jpl), INTENT(inout) :: pdevap_ice ! sublimation sensitivity
!
INTEGER :: jl ! dummy loop index
!
REAL(wp), DIMENSION(jpi,jpj) :: z1_at_i ! inverse of concentration
REAL(wp), DIMENSION(A2D(0)) :: z1_at_i ! inverse of concentration
!
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_qsr_m ! Mean solar heat flux over all categories
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_qns_m ! Mean non solar heat flux over all categories
......@@ -195,7 +235,7 @@ CONTAINS
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztem_m ! Mean temperature over all categories
!!----------------------------------------------------------------------
!
WHERE ( at_i (:,:) > 0._wp ) ; z1_at_i(:,:) = 1._wp / at_i (:,:)
WHERE ( pat_i(:,:) > 0._wp ) ; z1_at_i(:,:) = 1._wp / pat_i(:,:)
ELSEWHERE ; z1_at_i(:,:) = 0._wp
END WHERE
......@@ -203,13 +243,13 @@ CONTAINS
!
CASE( 0 , 1 )
!
ALLOCATE( z_qns_m(jpi,jpj), z_qsr_m(jpi,jpj), z_dqn_m(jpi,jpj), z_evap_m(jpi,jpj), z_devap_m(jpi,jpj) )
ALLOCATE( z_qns_m(A2D(0)), z_qsr_m(A2D(0)), z_dqn_m(A2D(0)), z_evap_m(A2D(0)), z_devap_m(A2D(0)) )
!
z_qns_m (:,:) = SUM( a_i(:,:,:) * pqns_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_qsr_m (:,:) = SUM( a_i(:,:,:) * pqsr_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_dqn_m (:,:) = SUM( a_i(:,:,:) * pdqn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
z_qns_m (:,:) = SUM( pa_i(:,:,:) * pqns_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_qsr_m (:,:) = SUM( pa_i(:,:,:) * pqsr_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_dqn_m (:,:) = SUM( pa_i(:,:,:) * pdqn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_evap_m (:,:) = SUM( pa_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_devap_m(:,:) = SUM( pa_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
DO jl = 1, jpl
pqns_ice (:,:,jl) = z_qns_m (:,:)
pqsr_ice (:,:,jl) = z_qsr_m (:,:)
......@@ -226,10 +266,10 @@ CONTAINS
!
CASE( 1 , 2 )
!
ALLOCATE( zalb_m(jpi,jpj), ztem_m(jpi,jpj) )
ALLOCATE( zalb_m(A2D(0)), ztem_m(A2D(0)) )
!
zalb_m(:,:) = SUM( a_i(:,:,:) * palb_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
ztem_m(:,:) = SUM( a_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
zalb_m(:,:) = SUM( pa_i(:,:,:) * palb_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
ztem_m(:,:) = SUM( pa_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
DO jl = 1, jpl
pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
......@@ -257,7 +297,7 @@ CONTAINS
REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos, zu_io, zv_io, zu_iom1, zv_iom1
REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04)
REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient
REAL(wp), DIMENSION(jpi,jpj) :: zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
REAL(wp), DIMENSION(A2D(0)) :: zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
!!-----------------------------------------------------------------------
!
! computation of friction velocity at T points
......@@ -268,36 +308,33 @@ CONTAINS
zv_io = v_ice(ji ,jj ) - ssv_m(ji ,jj )
zv_iom1 = v_ice(ji ,jj-1) - ssv_m(ji ,jj-1)
!
zfric(ji,jj) = rn_cio * ( 0.5_wp * ( zu_io*zu_io + zu_iom1*zu_iom1 + zv_io*zv_io + zv_iom1*zv_iom1 ) ) * tmask(ji,jj,1)
zfric(ji,jj) = rn_cio * ( 0.5_wp * ( zu_io*zu_io + zu_iom1*zu_iom1 + zv_io*zv_io + zv_iom1*zv_iom1 ) ) * smask0(ji,jj)
zvel (ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) + &
& ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )
END_2D
ELSE ! if no ice dynamics => transfer directly the atmospheric stress to the ocean
DO_2D( 0, 0, 0, 0 )
zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp * &
& ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) &
& + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
zvel(ji,jj) = 0._wp
zfric(ji,jj) = r1_rho0 * SQRT( utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj) ) * smask0(ji,jj)
zvel (ji,jj) = 0._wp
END_2D
ENDIF
CALL lbc_lnk( 'icesbc', zfric, 'T', 1.0_wp, zvel, 'T', 1.0_wp )
!
!--------------------------------------------------------------------!
! Partial computation of forcing for the thermodynamic sea ice model
!--------------------------------------------------------------------!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! needed for qlead
rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
DO_2D( 0, 0, 0, 0 ) ! needed for qlead
rswitch = smask0(ji,jj) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
!
! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- !
zqld = tmask(ji,jj,1) * rDt_ice * &
zqld = smask0(ji,jj) * rDt_ice * &
& ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) + &
& ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- !
! (mostly<0 but >0 if supercooling)
zqfr = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1) ! both < 0 (t_bo < sst) and > 0 (t_bo > sst)
zqfr_neg = MIN( zqfr , 0._wp ) ! only < 0
zqfr_pos = MAX( zqfr , 0._wp ) ! only > 0
zqfr = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * smask0(ji,jj) ! both < 0 (t_bo < sst) and > 0 (t_bo > sst)
zqfr_neg = MIN( zqfr , 0._wp ) ! only < 0
zqfr_pos = MAX( zqfr , 0._wp ) ! only > 0
! --- Sensible ocean-to-ice heat flux (W/m2) --- !
! (mostly>0 but <0 if supercooling)
......
src/ICE/icestp.F90
View file @ 3efbcb17
......@@ -77,6 +77,7 @@ MODULE icestp
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE timing ! Timing
USE prtctl ! Print control
USE lbclnk ! lateral boundary conditions (or mpp links)
IMPLICIT NONE
PRIVATE
......@@ -109,7 +110,7 @@ CONTAINS
!! - save the outputs
!! - save the outputs for restart when necessary
!!
!! ** Action : - time evolution of the LIM sea-ice model
!! ** Action : - time evolution of the SI3 sea-ice model
!! - update all sbc variables below sea-ice:
!! utau, vtau, taum, wndm, qns , qsr, emp , sfx
!!---------------------------------------------------------------------
......@@ -117,7 +118,7 @@ CONTAINS
INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices
INTEGER, INTENT(in) :: ksbc ! flux formulation (user defined, bulk, or Pure Coupled)
!
INTEGER :: jl ! dummy loop index
INTEGER :: ji, jj, jl ! dummy loop index
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('icestp')
......@@ -128,9 +129,12 @@ CONTAINS
!
kt_ice = kt ! -- Ice model time step
!
u_oce(:,:) = ssu_m(:,:) ! -- mean surface ocean current
v_oce(:,:) = ssv_m(:,:)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! -- mean surface ocean current
u_oce(ji,jj) = ssu_m(ji,jj)
v_oce(ji,jj) = ssv_m(ji,jj)
END_2D
!
! clem: I think t_bo needs to be defined everywhere but check
CALL eos_fzp( sss_m(:,:) , t_bo(:,:) ) ! -- freezing temperature [Kelvin] (set to rt0 over land)
t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
!
......@@ -152,6 +156,7 @@ CONTAINS
! utau_ice, vtau_ice = surface ice stress [N/m2]
!------------------------------------------------!
CALL ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice )
! => clem: here utau_ice and vtau_ice are defined everywhere
!-------------------------------------!
! --- ice dynamics and advection --- !
!-------------------------------------!
......@@ -161,6 +166,11 @@ CONTAINS
IF( ln_icedyn .AND. .NOT.ln_c1d ) &
& CALL ice_dyn( kt, Kmm ) ! -- Ice dynamics
!
! ==> clem: here, all the global variables are correctly defined in the halos:
! ato_i, a_i, v_i, v_s, sv_i, oa_i, a_ip, v_ip, v_il, e_i, e_s
! as well as h_i and t_su for practical reasons
! It may not be necessary. If it is not, then remove lbc_lnk in icedyn_rdgrft and iceitd_reb
CALL diag_trends( 1 ) ! record dyn trends
!
! !== lateral boundary conditions ==!
......@@ -170,7 +180,8 @@ CONTAINS
CALL ice_var_glo2eqv ! h_i and h_s for ice albedo calculation
CALL ice_var_agg(1) ! at_i for coupling
CALL store_fields ! Store now ice values
!
! ==> clem: here, full arrays = vt_i, vt_s, vt_ip, vt_il, at_i, at_ip, ato_i
! reduced arrays = st_i, et_i, et_s, tm_su
!------------------------------------------------------!
! --- Thermodynamical coupling with the atmosphere --- !
!------------------------------------------------------!
......@@ -184,11 +195,15 @@ CONTAINS
! qemp_oce, qemp_ice, = sensible heat (associated with evap & precip) [W/m2]
! qprec_ice, qevap_ice
!------------------------------------------------------!
! ==> clem: From here on, we only need to work on the interior domain
! though it necessitates a large lbc at the end of the time step
CALL ice_sbc_flx( kt, ksbc )
!----------------------------!
! --- ice thermodynamics --- !
!----------------------------!
IF( ln_icethd ) CALL ice_thd( kt ) ! -- Ice thermodynamics
!
! ==> clem: here, all the global variables are correctly defined in the halos
!
CALL diag_trends( 2 ) ! record thermo trends
CALL ice_var_glo2eqv ! necessary calls (at least for coupling)
......@@ -288,7 +303,7 @@ CONTAINS
CALL ice_drift_init ! initialization for diags of conservation
!
fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction
tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu
tn_ice(:,:,:) = t_su(A2D(0),:) ! initialisation of surface temp for coupled simu
!
IF( ln_rstart ) THEN
CALL iom_close( numrir ) ! close input ice restart file
......@@ -369,26 +384,33 @@ CONTAINS
INTEGER :: ji, jj, jl ! dummy loop index
!!----------------------------------------------------------------------
!
a_i_b (:,:,:) = a_i (:,:,:) ! ice area
v_i_b (:,:,:) = v_i (:,:,:) ! ice volume
v_s_b (:,:,:) = v_s (:,:,:) ! snow volume
v_ip_b(:,:,:) = v_ip(:,:,:) ! pond volume
v_il_b(:,:,:) = v_il(:,:,:) ! pond lid volume
sv_i_b(:,:,:) = sv_i(:,:,:) ! salt content
e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy
e_i_b (:,:,:,:) = e_i (:,:,:,:) ! ice thermal energy
WHERE( a_i_b(:,:,:) >= epsi20 )
h_i_b(:,:,:) = v_i_b(:,:,:) / a_i_b(:,:,:) ! ice thickness
h_s_b(:,:,:) = v_s_b(:,:,:) / a_i_b(:,:,:) ! snw thickness
ELSEWHERE
h_i_b(:,:,:) = 0._wp
h_s_b(:,:,:) = 0._wp
END WHERE
!
! ice velocities & total concentration
DO jl = 1, jpl
DO_2D( 0, 0, 0, 0 )
a_i_b (ji,jj,jl) = a_i (ji,jj,jl) ! ice area
v_i_b (ji,jj,jl) = v_i (ji,jj,jl) ! ice volume
v_s_b (ji,jj,jl) = v_s (ji,jj,jl) ! snow volume
v_ip_b(ji,jj,jl) = v_ip(ji,jj,jl) ! pond volume
v_il_b(ji,jj,jl) = v_il(ji,jj,jl) ! pond lid volume
sv_i_b(ji,jj,jl) = sv_i(ji,jj,jl) ! salt content
IF( a_i_b(ji,jj,jl) >= epsi20 ) THEN
h_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / a_i_b(ji,jj,jl) ! ice thickness
h_s_b(ji,jj,jl) = v_s_b(ji,jj,jl) / a_i_b(ji,jj,jl) ! snw thickness
ELSE
h_i_b(ji,jj,jl) = 0._wp
h_s_b(ji,jj,jl) = 0._wp
ENDIF
e_s_b (ji,jj,:,jl) = e_s (ji,jj,:,jl) ! snow thermal energy
e_i_b (ji,jj,:,jl) = e_i (ji,jj,:,jl) ! ice thermal energy
END_2D
ENDDO
! total concentration
at_i_b(:,:) = SUM( a_i_b(:,:,:), dim=3 )
u_ice_b(:,:) = u_ice(:,:)
v_ice_b(:,:) = v_ice(:,:)
! ice velocity
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
u_ice_b(ji,jj) = u_ice(ji,jj)
v_ice_b(ji,jj) = v_ice(ji,jj)
END_2D
!
END SUBROUTINE store_fields
......@@ -402,20 +424,23 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jl ! dummy loop index
!!----------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
sfx_res(ji,jj) = 0._wp ; wfx_res(ji,jj) = 0._wp ; hfx_res(ji,jj) = 0._wp
END_2D
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! needed for (at least) diag_adv_mass -> to be removed
DO_2D( 0, 0, 0, 0 )
sfx (ji,jj) = 0._wp ;
sfx_bri(ji,jj) = 0._wp ; sfx_lam(ji,jj) = 0._wp
sfx_sni(ji,jj) = 0._wp ; sfx_opw(ji,jj) = 0._wp
sfx_bog(ji,jj) = 0._wp ; sfx_dyn(ji,jj) = 0._wp
sfx_bom(ji,jj) = 0._wp ; sfx_sum(ji,jj) = 0._wp
sfx_res(ji,jj) = 0._wp ; sfx_sub(ji,jj) = 0._wp
sfx_sub(ji,jj) = 0._wp
!
wfx_snw(ji,jj) = 0._wp ; wfx_ice(ji,jj) = 0._wp
wfx_sni(ji,jj) = 0._wp ; wfx_opw(ji,jj) = 0._wp
wfx_bog(ji,jj) = 0._wp ; wfx_dyn(ji,jj) = 0._wp
wfx_bom(ji,jj) = 0._wp ; wfx_sum(ji,jj) = 0._wp
wfx_res(ji,jj) = 0._wp ; wfx_sub(ji,jj) = 0._wp
wfx_sub(ji,jj) = 0._wp
wfx_spr(ji,jj) = 0._wp ; wfx_lam(ji,jj) = 0._wp
wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp
wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp
......@@ -426,7 +451,7 @@ CONTAINS
hfx_snw(ji,jj) = 0._wp ; hfx_opw(ji,jj) = 0._wp
hfx_bog(ji,jj) = 0._wp ; hfx_dyn(ji,jj) = 0._wp
hfx_bom(ji,jj) = 0._wp ; hfx_sum(ji,jj) = 0._wp
hfx_res(ji,jj) = 0._wp ; hfx_sub(ji,jj) = 0._wp
hfx_sub(ji,jj) = 0._wp
hfx_spr(ji,jj) = 0._wp ; hfx_dif(ji,jj) = 0._wp
hfx_err_dif(ji,jj) = 0._wp
wfx_err_sub(ji,jj) = 0._wp
......@@ -451,7 +476,7 @@ CONTAINS
END_2D
DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
! SIMIP diagnostics
t_si (ji,jj,jl) = rt0 ! temp at the ice-snow interface
qcn_ice_bot(ji,jj,jl) = 0._wp
......@@ -477,22 +502,25 @@ CONTAINS
!! and outputs
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kn ! 1 = after dyn ; 2 = after thermo
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jl ! dummy loop index
!!----------------------------------------------------------------------
!
! --- trends of heat, salt, mass (used for conservation controls)
IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN
!
diag_heat(:,:) = diag_heat(:,:) &
& - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice &
& - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice
diag_sice(:,:) = diag_sice(:,:) &
& + SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice * rhoi
diag_vice(:,:) = diag_vice(:,:) &
& + SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhoi
diag_vsnw(:,:) = diag_vsnw(:,:) &
& + SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos
diag_vpnd(:,:) = diag_vpnd(:,:) &
& + SUM( v_ip + v_il - v_ip_b - v_il_b , dim=3 ) * r1_Dt_ice * rhow
DO_2D( 0, 0, 0, 0 )
diag_heat(ji,jj) = diag_heat(ji,jj) &
& - SUM(SUM( e_i (ji,jj,1:nlay_i,:) - e_i_b (ji,jj,1:nlay_i,:), dim=2 ) ) * r1_Dt_ice &
& - SUM(SUM( e_s (ji,jj,1:nlay_s,:) - e_s_b (ji,jj,1:nlay_s,:), dim=2 ) ) * r1_Dt_ice
diag_sice(ji,jj) = diag_sice(ji,jj) &
& + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * r1_Dt_ice * rhoi
diag_vice(ji,jj) = diag_vice(ji,jj) &
& + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * r1_Dt_ice * rhoi
diag_vsnw(ji,jj) = diag_vsnw(ji,jj) &
& + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * r1_Dt_ice * rhos
diag_vpnd(ji,jj) = diag_vpnd(ji,jj) &
& + SUM( v_ip(ji,jj,:)+v_il(ji,jj,:) - v_ip_b(ji,jj,:)-v_il_b(ji,jj,:) ) * r1_Dt_ice * rhow
END_2D
!
IF( kn == 2 ) CALL iom_put ( 'hfxdhc' , diag_heat ) ! output of heat trend
!
......@@ -501,10 +529,12 @@ CONTAINS
! --- trends of concentration (used for simip outputs)
IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN
!
diag_aice(:,:) = diag_aice(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice
DO_2D( 0, 0, 0, 0 )
diag_aice(ji,jj) = diag_aice(ji,jj) + SUM( a_i(ji,jj,:) - a_i_b(ji,jj,:) ) * r1_Dt_ice
END_2D
!
IF( kn == 1 ) CALL iom_put( 'afxdyn' , diag_aice ) ! dyn trend
IF( kn == 2 ) CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice ) ! thermo trend
IF( kn == 2 ) CALL iom_put( 'afxthd' , SUM( a_i(A2D(0),:) - a_i_b(A2D(0),:), dim=3 ) * r1_Dt_ice ) ! thermo trend
IF( kn == 2 ) CALL iom_put( 'afxtot' , diag_aice ) ! total trend
!
ENDIF
......
src/ICE/icetab.F90
View file @ 3efbcb17
......@@ -38,15 +38,24 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ndim1d ! 1d size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in ) :: tab2d ! input 2D field
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: tab2d ! input 2D field
REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) :: tab1d ! output 1D field
!
INTEGER :: jl, jn, jid, jjd
INTEGER :: ipi, ipj, ji0, jj0, jl, jn, jid, jjd
!!----------------------------------------------------------------------
ipi = SIZE(tab2d,1) ! 1st dimension
ipj = SIZE(tab2d,2) ! 2nd dimension
!
IF( ipi == jpi .AND. ipj == jpj ) THEN ! full arrays then no need to change index jid and jjd
ji0 = 0 ; jj0 = 0
ELSE ! reduced arrays then need to shift index by nn_hls
ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls
ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls)
!
DO jl = 1, jpl
DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1
jjd = ( tab_ind(jn) - 1 ) / jpi + 1
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
tab1d(jn,jl) = tab2d(jid,jjd,jl)
END DO
END DO
......@@ -59,14 +68,23 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ndim1d ! 1d size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: tab2d ! input 2D field
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: tab2d ! input 2D field
REAL(wp), DIMENSION(ndim1d) , INTENT(inout) :: tab1d ! output 1D field
!
INTEGER :: jn , jid, jjd
INTEGER :: ipi, ipj, ji0, jj0, jn, jid, jjd
!!----------------------------------------------------------------------
ipi = SIZE(tab2d,1) ! 1st dimension
ipj = SIZE(tab2d,2) ! 2nd dimension
!
IF( ipi == jpi .AND. ipj == jpj ) THEN ! full arrays then no need to change index jid and jjd
ji0 = 0 ; jj0 = 0
ELSE ! reduced arrays then need to shift index by nn_hls
ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls
ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls)
!
DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1
jjd = ( tab_ind(jn) - 1 ) / jpi + 1
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
tab1d( jn) = tab2d( jid, jjd)
END DO
END SUBROUTINE tab_2d_1d
......@@ -79,14 +97,23 @@ CONTAINS
INTEGER , INTENT(in ) :: ndim1d ! 1D size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in ) :: tab1d ! input 1D field
REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) :: tab2d ! output 2D field
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: tab2d ! output 2D field
!
INTEGER :: jl, jn, jid, jjd
INTEGER :: ipi, ipj, ji0, jj0, jl, jn, jid, jjd
!!----------------------------------------------------------------------
ipi = SIZE(tab2d,1) ! 1st dimension
ipj = SIZE(tab2d,2) ! 2nd dimension
!
IF( ipi == jpi .AND. ipj == jpj ) THEN ! full arrays then no need to change index jid and jjd
ji0 = 0 ; jj0 = 0
ELSE ! reduced arrays then need to shift index by nn_hls
ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls
ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls)
!
DO jl = 1, jpl
DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1
jjd = ( tab_ind(jn) - 1 ) / jpi + 1
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
tab2d(jid,jjd,jl) = tab1d(jn,jl)
END DO
END DO
......@@ -100,13 +127,22 @@ CONTAINS
INTEGER , INTENT(in ) :: ndim1d ! 1D size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
REAL(wp), DIMENSION(ndim1d) , INTENT(in ) :: tab1d ! input 1D field
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: tab2d ! output 2D field
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: tab2d ! output 2D field
!
INTEGER :: jn , jid, jjd
INTEGER :: ipi, ipj, ji0, jj0, jn, jid, jjd
!!----------------------------------------------------------------------
ipi = SIZE(tab2d,1) ! 1st dimension
ipj = SIZE(tab2d,2) ! 2nd dimension
!
IF( ipi == jpi .AND. ipj == jpj ) THEN ! full arrays then no need to change index jid and jjd
ji0 = 0 ; jj0 = 0
ELSE ! reduced arrays then need to shift index by nn_hls
ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls
ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls)
!
DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1
jjd = ( tab_ind(jn) - 1 ) / jpi + 1
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
tab2d(jid, jjd) = tab1d( jn)
END DO
END SUBROUTINE tab_1d_2d
......
src/ICE/icethd.F90
View file @ 3efbcb17
......@@ -98,7 +98,7 @@ CONTAINS
! convergence tests
IF( ln_zdf_chkcvg ) THEN
ALLOCATE( ztice_cvgerr(jpi,jpj,jpl) , ztice_cvgstp(jpi,jpj,jpl) )
ALLOCATE( ztice_cvgerr(A2D(0),jpl) , ztice_cvgstp(A2D(0),jpl) )
ztice_cvgerr = 0._wp ; ztice_cvgstp = 0._wp
ENDIF
!
......@@ -111,7 +111,7 @@ CONTAINS
! select ice covered grid points
npti = 0 ; nptidx(:) = 0
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF ( a_i(ji,jj,jl) > epsi10 ) THEN
npti = npti + 1
nptidx(npti) = (jj - 1) * jpi + ji
......@@ -158,6 +158,13 @@ CONTAINS
IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- !
!
IF( ln_icedO ) CALL ice_thd_do ! --- Frazil ice growth in leads --- !
!
! ! --- LBC for the halos --- !
! the 2 lbc below could be avoided if calculations above were performed over the full domain
! but we think it is more efficient this way
CALL lbc_lnk( 'icethd', a_i , 'T', 1._wp, v_i , 'T', 1._wp, v_s , 'T', 1._wp, sv_i, 'T', 1._wp, oa_i, 'T', 1._wp, &
& a_ip, 'T', 1._wp, v_ip, 'T', 1._wp, v_il, 'T', 1._wp, t_su, 'T', 1._wp )
CALL lbc_lnk( 'icethd', e_i , 'T', 1._wp, e_s , 'T', 1._wp )
!
CALL ice_cor( kt , 2 ) ! --- Corrections --- !
!
......
src/ICE/icethd_do.F90
View file @ 3efbcb17
......@@ -28,7 +28,6 @@ MODULE icethd_do
USE in_out_manager ! I/O manager
USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE lbclnk ! lateral boundary conditions (or mpp links)
IMPLICIT NONE
PRIVATE
......@@ -105,7 +104,7 @@ CONTAINS
IF( ln_icediachk ) CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
IF( ln_icediachk ) CALL ice_cons2D ( 0, 'icethd_do', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft )
at_i(:,:) = SUM( a_i, dim=3 )
at_i(A2D(0)) = SUM( a_i(A2D(0),:), dim=3 )
!------------------------------------------------------------------------------!
! 1) Compute thickness, salinity, enthalpy, age, area and volume of new ice
!------------------------------------------------------------------------------!
......@@ -113,7 +112,7 @@ CONTAINS
! Identify grid points where new ice forms
npti = 0 ; nptidx(:) = 0
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF ( qlead(ji,jj) < 0._wp ) THEN
npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji
......@@ -325,6 +324,8 @@ CONTAINS
!
ENDIF ! npti > 0
!
! the following fields need to be updated on the halos (done in icethd): a_i, v_i, sv_i, e_i
!
IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
IF( ln_icediachk ) CALL ice_cons2D (1, 'icethd_do', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft)
!
......@@ -372,8 +373,8 @@ CONTAINS
DO_2D( 0, 0, 0, 0 )
IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling
! -- Wind stress -- !
ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) + utau_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + vtau_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
ztaux = utau_ice(ji,jj) * smask0(ji,jj)
ztauy = vtau_ice(ji,jj) * smask0(ji,jj)
! Square root of wind stress
ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
......@@ -415,8 +416,6 @@ CONTAINS
!
END_2D
!
CALL lbc_lnk( 'icethd_frazil', fraz_frac, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp )
ENDIF
END SUBROUTINE ice_thd_frazil
......