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
Expand all 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 ...@@ -44,6 +44,8 @@ MODULE sbcabl
PUBLIC sbc_abl_init ! routine called in sbcmod module PUBLIC sbc_abl_init ! routine called in sbcmod module
PUBLIC sbc_abl ! 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) !! NEMO/OPA 3.7 , NEMO-consortium (2014)
!! $Id: sbcabl.F90 6416 2016-04-01 12:22:17Z clem $ !! $Id: sbcabl.F90 6416 2016-04-01 12:22:17Z clem $
...@@ -308,8 +310,8 @@ CONTAINS ...@@ -308,8 +310,8 @@ CONTAINS
!! - Perform 1 time-step of the ABL model !! - Perform 1 time-step of the ABL model
!! - Finalize flux computation in blk_oce_2 !! - Finalize flux computation in blk_oce_2
!! !!
!! ** Outputs : - utau : i-component of the stress at U-point (N/m2) !! ** Outputs : - utau : i-component of the stress at T-point (N/m2)
!! - vtau : j-component of the stress at V-point (N/m2) !! - vtau : j-component of the stress at T-point (N/m2)
!! - taum : Wind stress module at T-point (N/m2) !! - taum : Wind stress module at T-point (N/m2)
!! - wndm : Wind speed module at T-point (m/s) !! - wndm : Wind speed module at T-point (m/s)
!! - qsr : Solar heat flux over the ocean (W/m2) !! - qsr : Solar heat flux over the ocean (W/m2)
...@@ -340,7 +342,7 @@ CONTAINS ...@@ -340,7 +342,7 @@ CONTAINS
CALL blk_oce_1( kt, u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in 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 & 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_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in
& sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) , & ! <<= in & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) , & ! <<= in
& tsk_m, zssq, zcd_du, zsen, zlat, zevp ) ! =>> out & tsk_m, zssq, zcd_du, zsen, zlat, zevp ) ! =>> out
...@@ -348,7 +350,7 @@ CONTAINS ...@@ -348,7 +350,7 @@ CONTAINS
#if defined key_si3 #if defined key_si3
CALL blk_ice_1( u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in 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 & 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 & pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui ) ! <<= out
#endif #endif
......
src/ICE/ice.F90
View file @ 3efbcb17
...@@ -451,6 +451,8 @@ MODULE ice ...@@ -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_bot !: Bottom conduction flux (W/m2)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qcn_ice_top !: Surface 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) !! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: ice.F90 15388 2021-10-17 11:33:47Z clem $ !! $Id: ice.F90 15388 2021-10-17 11:33:47Z clem $
...@@ -464,71 +466,103 @@ CONTAINS ...@@ -464,71 +466,103 @@ CONTAINS
!!----------------------------------------------------------------- !!-----------------------------------------------------------------
INTEGER :: ice_alloc INTEGER :: ice_alloc
! !
INTEGER :: ierr(16), ii INTEGER :: ierr(21), ii
!!----------------------------------------------------------------- !!-----------------------------------------------------------------
ierr(:) = 0 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 ii = ii + 1
ALLOCATE( u_oce (jpi,jpj) , v_oce (jpi,jpj) , ht_i_new (jpi,jpj) , fraz_frac (jpi,jpj) , & ALLOCATE( h_i (jpi,jpj,jpl) , a_i (jpi,jpj,jpl) , v_i (jpi,jpj,jpl) , &
& strength (jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) , & & v_s (jpi,jpj,jpl) , h_s (jpi,jpj,jpl) , &
& delta_i (jpi,jpj) , divu_i (jpi,jpj) , shear_i (jpi,jpj) , & & s_i (jpi,jpj,jpl) , sv_i(jpi,jpj,jpl) , o_i (jpi,jpj,jpl) , oa_i (jpi,jpj,jpl) , &
& aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv (jpi,jpj) , STAT=ierr(ii) ) & 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 ii = ii + 1
ALLOCATE( t_bo (jpi,jpj) , wfx_snw_sni(jpi,jpj) , & ALLOCATE( e_s(jpi,jpj,nlay_s,jpl) , e_i(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) )
& 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) )
! * Ice global state variables ! * Before values of global variables
ii = ii + 1 ii = ii + 1
ALLOCATE( qtr_ice_bot(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) , & ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) )
& 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) )
! * fluxes
ii = ii + 1 ii = ii + 1
ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , & 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
& 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) , & ! * ice rheology
& sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s(jpi,jpj) , & ii = ii+1
& om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj), icb_mask(jpi,jpj), STAT=ierr(ii) ) 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 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 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 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), & 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_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , & & v_s_b (A2D(0),jpl) , h_s_b (A2D(0),jpl) , &
& dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) , STAT = ierr(ii) ) & 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 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 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), & ALLOCATE( st_i (A2D(0)) , et_i (A2D(0)) , et_s(A2D(0)) , hm_i (A2D(0)) , &
& v_ip_b(jpi,jpj,jpl) , v_il_b(jpi,jpj,jpl) , & & hm_ip(A2D(0)) , hm_il(A2D(0)) , tm_i(A2D(0)) , tm_s (A2D(0)) , &
& 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) , & & sm_i (A2D(0)) , hm_s (A2D(0)) , om_i(A2D(0)) , bvm_i(A2D(0)) , &
& STAT=ierr(ii) ) & tm_su(A2D(0)) , STAT=ierr(ii) )
! * others
ii = ii + 1 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 ! * Ice thickness distribution variables
ii = ii + 1 ii = ii + 1
...@@ -536,19 +570,19 @@ CONTAINS ...@@ -536,19 +570,19 @@ CONTAINS
! * Ice diagnostics ! * Ice diagnostics
ii = ii + 1 ii = ii + 1
ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & ALLOCATE( diag_trp_vi (A2D(0)) , diag_trp_vs (A2D(0)) , diag_trp_ei (A2D(0)) , &
& diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), & & diag_trp_es (A2D(0)) , diag_trp_sv (A2D(0)) , diag_heat (A2D(0)) , &
& diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), diag_vpnd(jpi,jpj), & & diag_sice (A2D(0)) , diag_vice (A2D(0)) , diag_vsnw (A2D(0)) , diag_aice(A2D(0)) , diag_vpnd(A2D(0)), &
& diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) & diag_adv_mass(A2D(0)) , diag_adv_salt(A2D(0)) , diag_adv_heat(A2D(0)) , STAT=ierr(ii) )
! * Ice conservation ! * Ice conservation
ii = ii + 1 ii = ii + 1
ALLOCATE( diag_v (jpi,jpj) , diag_s (jpi,jpj) , diag_t (jpi,jpj), & ALLOCATE( diag_v (A2D(0)) , diag_s (A2D(0)) , diag_t (A2D(0)), &
& diag_fv(jpi,jpj) , diag_fs(jpi,jpj) , diag_ft(jpi,jpj), STAT=ierr(ii) ) & diag_fv(A2D(0)) , diag_fs(A2D(0)) , diag_ft(A2D(0)), STAT=ierr(ii) )
! * SIMIP diagnostics ! * SIMIP diagnostics
ii = ii + 1 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(:) ) ice_alloc = MAXVAL( ierr(:) )
IF( ice_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_alloc: failed to allocate arrays.' ) 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 ...@@ -48,7 +48,7 @@ MODULE icealb
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
CONTAINS 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 *** !! *** ROUTINE ice_alb ***
!! !!
...@@ -94,16 +94,16 @@ CONTAINS ...@@ -94,16 +94,16 @@ CONTAINS
!! Brandt et al. 2005, J. Climate, vol 18 !! Brandt et al. 2005, J. Climate, vol 18
!! Grenfell & Perovich 2004, JGR, vol 109 !! Grenfell & Perovich 2004, JGR, vol 109
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: pt_su ! ice surface temperature (Kelvin) LOGICAL , INTENT(in ) :: ld_pnd_alb ! effect of melt ponds on albedo
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: ph_ice ! sea-ice thickness REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: pt_su ! ice surface temperature (Kelvin)
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: ph_snw ! snow depth REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: ph_ice ! sea-ice thickness
LOGICAL , INTENT(in ) :: ld_pnd_alb ! effect of melt ponds on albedo REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: ph_snw ! snow depth
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: pafrac_pnd ! melt pond relative fraction (per unit ice area) REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: 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(A2D(0),jpl) :: ph_pnd ! melt pond depth
REAL(wp), INTENT(in ), DIMENSION(:,:) :: pcloud_fra ! cloud fraction REAL(wp), INTENT(in ), DIMENSION(A2D(0)) :: pcloud_fra ! cloud fraction
REAL(wp), INTENT( out), DIMENSION(:,:,:) :: palb_ice ! albedo of ice 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 INTEGER :: ji, jj, jl ! dummy loop indices
REAL(wp) :: z1_c1, z1_c2,z1_c3, z1_c4 ! local scalar 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) REAL(wp) :: z1_href_pnd ! inverse of the characteristic length scale (Lecomte et al. 2015)
...@@ -121,10 +121,10 @@ CONTAINS ...@@ -121,10 +121,10 @@ CONTAINS
z1_c3 = 1._wp / 0.02_wp z1_c3 = 1._wp / 0.02_wp
z1_c4 = 1._wp / 0.03_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 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 ---! !--- Specific snow, ice and pond fractions ---!
...@@ -164,10 +164,10 @@ CONTAINS ...@@ -164,10 +164,10 @@ CONTAINS
zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd(ji,jj,jl) * z1_href_pnd ) 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 ! !--- 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 & 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 ! 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 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 ...@@ -83,25 +83,33 @@ CONTAINS
CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 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 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) :: zdiag_mass, zdiag_salt, zdiag_heat
REAL(wp), DIMENSION(jpi,jpj,10) :: ztmp3 REAL(wp), DIMENSION(A2D(0),10) :: ztmp3
REAL(wp), DIMENSION(jpi,jpj,jpl,8) :: ztmp4 REAL(wp), DIMENSION(A2D(0),jpl,8) :: ztmp4
REAL(wp), DIMENSION(10) :: zchk3 REAL(wp), DIMENSION(10) :: zchk3
REAL(wp), DIMENSION(8) :: zchk4 REAL(wp), DIMENSION(8) :: zchk4
!!------------------------------------------------------------------- !!-------------------------------------------------------------------
! !
! -- quantities -- ! DO_2D( 0, 0, 0, 0 )
ztmp3(:,:,1) = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ! volume ! -- quantities -- !
ztmp3(:,:,2) = SUM( sv_i * rhoi, dim=3 ) * e1e2t ! salt ztmp3(ji,jj,1) = SUM( v_i(ji,jj,:) * rhoi + v_s(ji,jj,:) * rhos + &
ztmp3(:,:,3) = ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ! heat & ( 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
! -- fluxes -- ! ztmp3(ji,jj,3) = ( SUM( SUM( e_i(ji,jj,:,:), dim=2 ) ) + & ! heat
ztmp3(:,:,4) = ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd & ! mass & SUM( SUM( e_s(ji,jj,:,:), dim=2 ) ) ) * e1e2t(ji,jj)
& + 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 ! -- fluxes -- !
& + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ztmp3(ji,jj,4) = ( wfx_bog (ji,jj) + wfx_bom (ji,jj) + wfx_sum (ji,jj) + wfx_sni (ji,jj) & ! mass
ztmp3(:,:,6) = ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & ! heat & + wfx_opw (ji,jj) + wfx_res (ji,jj) + wfx_dyn (ji,jj) + wfx_lam (ji,jj) + wfx_pnd(ji,jj) &
& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t & + 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 -- ! ! -- global sum -- !
zchk3(1:6) = glob_sum_vec( 'icectl', ztmp3(:,:,1:6) ) zchk3(1:6) = glob_sum_vec( 'icectl', ztmp3(:,:,1:6) )
...@@ -123,25 +131,33 @@ CONTAINS ...@@ -123,25 +131,33 @@ CONTAINS
zdiag_heat = ( zchk3(3) - pdiag_t ) * r1_Dt_ice + ( zchk3(6) - pdiag_ft ) zdiag_heat = ( zchk3(3) - pdiag_t ) * r1_Dt_ice + ( zchk3(6) - pdiag_ft )
! -- max concentration diag -- ! ! -- max concentration diag -- !
ztmp3(:,:,7) = SUM( a_i, dim=3 ) DO_2D( 0, 0, 0, 0 )
zchk3(7) = glob_max( 'icectl', ztmp3(:,:,7) ) ztmp3(ji,jj,7) = SUM( a_i(ji,jj,:) )
END_2D
zchk3(7) = glob_max( 'icectl', ztmp3(:,:,7) )
! -- advection scheme is conservative? -- ! ! -- advection scheme is conservative? -- !
ztmp3(:,:,8 ) = diag_adv_mass * e1e2t DO_2D( 0, 0, 0, 0 )
ztmp3(:,:,9 ) = diag_adv_heat * e1e2t ztmp3(ji,jj,8 ) = diag_adv_mass(ji,jj) * e1e2t(ji,jj)
ztmp3(:,:,10) = SUM( a_i + epsi10, dim=3 ) * e1e2t ! ice area (+epsi10 to set a threshold > 0 when there is no ice) ztmp3(ji,jj,9 ) = diag_adv_heat(ji,jj) * e1e2t(ji,jj)
zchk3(8:10) = glob_sum_vec( 'icectl', ztmp3(:,:,8:10) ) 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 -- ! ! -- min diags -- !
ztmp4(:,:,:,1) = v_i DO jl = 1, jpl
ztmp4(:,:,:,2) = v_s DO_2D( 0, 0, 0, 0 )
ztmp4(:,:,:,3) = v_ip ztmp4(ji,jj,jl,1) = v_i(ji,jj,jl)
ztmp4(:,:,:,4) = v_il ztmp4(ji,jj,jl,2) = v_s(ji,jj,jl)
ztmp4(:,:,:,5) = a_i ztmp4(ji,jj,jl,3) = v_ip(ji,jj,jl)
ztmp4(:,:,:,6) = sv_i ztmp4(ji,jj,jl,4) = v_il(ji,jj,jl)
ztmp4(:,:,:,7) = SUM( e_i, dim=3 ) ztmp4(ji,jj,jl,5) = a_i(ji,jj,jl)
ztmp4(:,:,:,8) = SUM( e_s, dim=3 ) ztmp4(ji,jj,jl,6) = sv_i(ji,jj,jl)
zchk4(1:8) = glob_min_vec( 'icectl', ztmp4(:,:,:,1:8) ) 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 IF( lwp ) THEN
! check conservation issues ! check conservation issues
...@@ -188,17 +204,21 @@ CONTAINS ...@@ -188,17 +204,21 @@ CONTAINS
!!------------------------------------------------------------------- !!-------------------------------------------------------------------
CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
!! !!
REAL(wp), DIMENSION(jpi,jpj,4) :: ztmp INTEGER :: ji, jj ! dummy loop index
REAL(wp), DIMENSION(4) :: zchk REAL(wp), DIMENSION(A2D(0),4) :: ztmp
REAL(wp), DIMENSION(4) :: zchk
!!------------------------------------------------------------------- !!-------------------------------------------------------------------
DO_2D( 0, 0, 0, 0 )
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(ji,jj,1) = ( wfx_ice (ji,jj) + wfx_snw (ji,jj) + wfx_pnd (ji,jj) + wfx_spr(ji,jj) + wfx_sub(ji,jj) &
ztmp(:,:,3) = ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ! heat & + diag_vice(ji,jj) + diag_vsnw(ji,jj) + diag_vpnd(ji,jj) - diag_adv_mass(ji,jj) ) * e1e2t(ji,jj) ! mass diag
! equivalent to this: ztmp(ji,jj,2) = ( sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) ) * e1e2t(ji,jj) ! salt
!! ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 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
!! & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t ) ! equivalent to this:
ztmp(:,:,4) = SUM( a_i + epsi10, dim=3 ) * e1e2t ! ice area (+epsi10 to set a threshold > 0 when there is no ice) !! ( -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 ! global sums
zchk(1:4) = glob_sum_vec( 'icectl', ztmp(:,:,1:4) ) zchk(1:4) = glob_sum_vec( 'icectl', ztmp(:,:,1:4) )
...@@ -226,11 +246,11 @@ CONTAINS ...@@ -226,11 +246,11 @@ CONTAINS
!!------------------------------------------------------------------- !!-------------------------------------------------------------------
INTEGER , INTENT(in) :: icount ! called at: =0 the begining of the routine, =1 the end 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 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, & REAL(wp), DIMENSION(A2D(0)) :: zdiag_mass, zdiag_salt, zdiag_heat, &
& zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax & zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax
INTEGER :: jl, jk INTEGER :: ji, jj, jl, jk
LOGICAL :: ll_stop_m = .FALSE. LOGICAL :: ll_stop_m = .FALSE.
LOGICAL :: ll_stop_s = .FALSE. LOGICAL :: ll_stop_s = .FALSE.
LOGICAL :: ll_stop_t = .FALSE. LOGICAL :: ll_stop_t = .FALSE.
...@@ -239,62 +259,79 @@ CONTAINS ...@@ -239,62 +259,79 @@ CONTAINS
! !
IF( icount == 0 ) THEN IF( icount == 0 ) THEN
pdiag_v = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) DO_2D( 0, 0, 0, 0 )
pdiag_s = SUM( sv_i * rhoi , dim=3 ) 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_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) 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 ! mass flux
pdiag_fv = wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & pdiag_fv(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) &
& wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr & + 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 ! 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 ! heat flux
pdiag_ft = hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 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 - hfx_dyn - hfx_res - hfx_sub - hfx_spr & - 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 ELSEIF( icount == 1 ) THEN
! -- mass diag -- ! ! -- mass diag -- !
zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) - pdiag_v ) * r1_Dt_ice & DO_2D( 0, 0, 0, 0 )
& + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & zdiag_mass(ji,jj) = ( SUM( v_i(ji,jj,:) * rhoi + v_s(ji,jj,:) * rhos &
& wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) & & + ( v_ip(ji,jj,:) + v_il(ji,jj,:) ) * rhow ) - pdiag_v(ji,jj) ) * r1_Dt_ice &
& - pdiag_fv & + ( 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. IF( MAXVAL( ABS(zdiag_mass) ) > rchk_m * rn_icechk_cel ) ll_stop_m = .TRUE.
! !
! -- salt diag -- ! ! -- salt diag -- !
zdiag_salt = ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_Dt_ice & DO_2D( 0, 0, 0, 0 )
& + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) & zdiag_salt(ji,jj) = ( SUM( sv_i(ji,jj,:) * rhoi ) - pdiag_s(ji,jj) ) * r1_Dt_ice &
& - pdiag_fs & + ( 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. IF( MAXVAL( ABS(zdiag_salt) ) > rchk_s * rn_icechk_cel ) ll_stop_s = .TRUE.
! !
! -- heat diag -- ! ! -- 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 & DO_2D( 0, 0, 0, 0 )
& + ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & zdiag_heat(ji,jj) = ( SUM( SUM( e_i(ji,jj,:,:), dim=2 ) ) &
& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) & & + SUM( SUM( e_s(ji,jj,:,:), dim=2 ) ) - pdiag_t(ji,jj) ) * r1_Dt_ice &
& - pdiag_ft & + ( 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. IF( MAXVAL( ABS(zdiag_heat) ) > rchk_t * rn_icechk_cel ) ll_stop_t = .TRUE.
! !
! -- other diags -- ! ! -- other diags -- !
! a_i < 0 ! a_i < 0
zdiag_amin(:,:) = 0._wp zdiag_amin(:,:) = 0._wp
DO jl = 1, jpl 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 ENDDO
! v_i < 0 ! v_i < 0
zdiag_vmin(:,:) = 0._wp zdiag_vmin(:,:) = 0._wp
DO jl = 1, jpl 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 ENDDO
! s_i < 0 ! s_i < 0
zdiag_smin(:,:) = 0._wp zdiag_smin(:,:) = 0._wp
DO jl = 1, jpl 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 ENDDO
! e_i < 0 ! e_i < 0
zdiag_emin(:,:) = 0._wp zdiag_emin(:,:) = 0._wp
DO jl = 1, jpl DO jl = 1, jpl
DO jk = 1, nlay_i 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
ENDDO ENDDO
! a_i > amax ! a_i > amax
...@@ -528,8 +565,8 @@ CONTAINS ...@@ -528,8 +565,8 @@ CONTAINS
INTEGER :: jl, ji, jj INTEGER :: jl, ji, jj
!!------------------------------------------------------------------- !!-------------------------------------------------------------------
DO ji = mi0(ki), mi1(ki) DO ji = mi0(ki,nn_hls), mi1(ki,nn_hls)
DO jj = mj0(kj), mj1(kj) DO jj = mj0(kj,nn_hls), mj1(kj,nn_hls)
WRITE(numout,*) ' time step ',kt,' ',cd1 ! print title WRITE(numout,*) ' time step ',kt,' ',cd1 ! print title
...@@ -733,10 +770,10 @@ CONTAINS ...@@ -733,10 +770,10 @@ CONTAINS
CALL prt_ctl_info(' ') CALL prt_ctl_info(' ')
CALL prt_ctl_info(' - Stresses : ') CALL prt_ctl_info(' - Stresses : ')
CALL prt_ctl_info(' ~~~~~~~~~~ ') CALL prt_ctl_info(' ~~~~~~~~~~ ')
CALL prt_ctl(tab2d_1=utau , clinfo1= ' utau : ', mask1 = umask, & CALL prt_ctl(tab2d_1=utau , clinfo1= ' utau : ', mask1 = tmask, &
& tab2d_2=vtau , clinfo2= ' vtau : ', mask2 = vmask) & tab2d_2=vtau , clinfo2= ' vtau : ', mask2 = tmask)
CALL prt_ctl(tab2d_1=utau_ice , clinfo1= ' utau_ice : ', mask1 = umask, & CALL prt_ctl(tab2d_1=utau_ice , clinfo1= ' utau_ice : ', mask1 = tmask, &
& tab2d_2=vtau_ice , clinfo2= ' vtau_ice : ', mask2 = vmask) & tab2d_2=vtau_ice , clinfo2= ' vtau_ice : ', mask2 = tmask)
END SUBROUTINE ice_prt3D END SUBROUTINE ice_prt3D
...@@ -751,8 +788,9 @@ CONTAINS ...@@ -751,8 +788,9 @@ CONTAINS
!!------------------------------------------------------------------- !!-------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ice time-step index INTEGER, INTENT(in) :: kt ! ice time-step index
! !
REAL(wp), DIMENSION(jpi,jpj,6) :: ztmp INTEGER :: ji, jj ! dummy loop index
REAL(wp), DIMENSION(6) :: zchk REAL(wp), DIMENSION(A2D(0),6) :: ztmp
REAL(wp), DIMENSION(6) :: zchk
!!------------------------------------------------------------------- !!-------------------------------------------------------------------
! !
IF( kt == nit000 .AND. lwp ) THEN IF( kt == nit000 .AND. lwp ) THEN
...@@ -762,25 +800,27 @@ CONTAINS ...@@ -762,25 +800,27 @@ CONTAINS
ENDIF ENDIF
! !
! -- 2D budgets (must be close to 0) -- ! ! -- 2D budgets (must be close to 0) -- !
ztmp(:,:,1) = wfx_ice (:,:) + wfx_snw (:,:) + wfx_spr (:,:) + wfx_sub(:,:) + wfx_pnd(:,:) & DO_2D( 0, 0, 0, 0 )
& + diag_vice(:,:) + diag_vsnw(:,:) + diag_vpnd(:,:) - diag_adv_mass(:,:) ztmp(ji,jj,1) = wfx_ice (ji,jj) + wfx_snw (ji,jj) + wfx_spr (ji,jj) + wfx_sub(ji,jj) + wfx_pnd(ji,jj) &
ztmp(:,:,2) = sfx(:,:) + diag_sice(:,:) - diag_adv_salt(:,:) & + diag_vice(ji,jj) + diag_vsnw(ji,jj) + diag_vpnd(ji,jj) - diag_adv_mass(ji,jj)
ztmp(:,:,3) = qt_oce_ai(:,:) - qt_atm_oi(:,:) + diag_heat(:,:) - diag_adv_heat(:,:) 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 ! write outputs
CALL iom_put( 'icedrift_mass', ztmp(:,:,1) ) CALL iom_put( 'icedrift_mass', ztmp(:,:,1) )
CALL iom_put( 'icedrift_salt', ztmp(:,:,2) ) CALL iom_put( 'icedrift_salt', ztmp(:,:,2) )
CALL iom_put( 'icedrift_heat', ztmp(:,:,3) ) CALL iom_put( 'icedrift_heat', ztmp(:,:,3) )
! -- 1D budgets -- ! ! -- 1D budgets -- !
ztmp(:,:,1) = ztmp(:,:,1) * e1e2t * rDt_ice ! mass DO_2D( 0, 0, 0, 0 )
ztmp(:,:,2) = ztmp(:,:,2) * e1e2t * rDt_ice * 1.e-3 ! salt ztmp(ji,jj,1) = ztmp(ji,jj,1) * e1e2t(ji,jj) * rDt_ice ! mass
ztmp(:,:,3) = ztmp(:,:,3) * e1e2t ! heat 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(:,:,4) = diag_adv_mass * e1e2t * rDt_ice
ztmp(:,:,5) = diag_adv_salt * e1e2t * rDt_ice * 1.e-3 ztmp(ji,jj,4) = diag_adv_mass(ji,jj) * e1e2t(ji,jj) * rDt_ice
ztmp(:,:,6) = diag_adv_heat * e1e2t 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 ! global sums
zchk(1:6) = glob_sum_vec( 'icectl', ztmp(:,:,1:6) ) zchk(1:6) = glob_sum_vec( 'icectl', ztmp(:,:,1:6) )
......
src/ICE/icedia.F90
View file @ 3efbcb17
...@@ -33,6 +33,9 @@ MODULE icedia ...@@ -33,6 +33,9 @@ MODULE icedia
PUBLIC ice_dia ! called by icestp.F90 PUBLIC ice_dia ! called by icestp.F90
PUBLIC ice_dia_init ! called in 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), 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), 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 REAL(wp) :: frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot ! global forcing trends
...@@ -48,7 +51,7 @@ CONTAINS ...@@ -48,7 +51,7 @@ CONTAINS
!!---------------------------------------------------------------------! !!---------------------------------------------------------------------!
!! *** ROUTINE ice_dia_alloc *** !! *** 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 ) CALL mpp_sum ( 'icedia', ice_dia_alloc )
IF( ice_dia_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_dia_alloc: failed to allocate arrays' ) IF( ice_dia_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_dia_alloc: failed to allocate arrays' )
...@@ -64,8 +67,9 @@ CONTAINS ...@@ -64,8 +67,9 @@ CONTAINS
!!--------------------------------------------------------------------------- !!---------------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time step INTEGER, INTENT(in) :: kt ! ocean time step
!! !!
REAL(wp), DIMENSION(jpi,jpj,16) :: ztmp INTEGER :: ji, jj ! dummy loop index
REAL(wp), DIMENSION(16) :: zbg REAL(wp), DIMENSION(A2D(0),16) :: ztmp
REAL(wp), DIMENSION(16) :: zbg
!!--------------------------------------------------------------------------- !!---------------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('ice_dia') IF( ln_timing ) CALL timing_start('ice_dia')
...@@ -85,31 +89,32 @@ CONTAINS ...@@ -85,31 +89,32 @@ CONTAINS
! 1 - Trends due to forcing ! ! 1 - Trends due to forcing !
! ---------------------------! ! ---------------------------!
! they must be kept outside an IF(iom_use) because of the call to dia_rst below ! 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 DO_2D( 0, 0, 0, 0 )
ztmp(:,:,2) = - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ! freshwater flux ice/snow-atm 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(:,:,3) = - sfx (:,:) * e1e2t(:,:) ! salt fluxes ice/snow-ocean ztmp(ji,jj,2) = - ( wfx_sub(ji,jj) + wfx_spr(ji,jj) ) * e1e2t(ji,jj) ! freshwater flux ice/snow-atm
ztmp(:,:,4) = qt_atm_oi(:,:) * e1e2t(:,:) ! heat on top of ice-ocean ztmp(ji,jj,3) = - sfx (ji,jj) * e1e2t(ji,jj) ! salt fluxes ice/snow-ocean
ztmp(:,:,5) = qt_oce_ai(:,:) * e1e2t(:,:) ! heat on top of ocean (and below ice) 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 ! ! 2 - Contents !
! ----------------------- ! ! ----------------------- !
IF( iom_use('ibgvol_tot' ) ) ztmp(:,:,6 ) = vt_i (:,:) * e1e2t(:,:) ! ice 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 (:,:) * e1e2t(:,:) ! snow 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 (:,:) * e1e2t(:,:) ! area 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(:,:) ! salt content IF( iom_use('ibgsalt_tot') ) ztmp(:,:,9 ) = st_i (:,:) * e1e2t(A2D(0)) ! salt content
IF( iom_use('ibgheat_tot') ) ztmp(:,:,10) = et_i (:,:) * e1e2t(:,:) ! heat 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(:,:) ! 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(:,:) * e1e2t(:,:) ! ice pond volume 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(:,:) * e1e2t(:,:) ! ice pond lid volume IF( iom_use('ilbgvol_tot') ) ztmp(:,:,13) = vt_il(A2D(0)) * e1e2t(A2D(0)) ! ice pond lid volume
! ---------------------------------- ! ! ---------------------------------- !
! 3 - Content variations and drifts ! ! 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('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(:,:) ! salt content 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') ) & 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 ! global sum
zbg(1:16) = glob_sum_vec( 'icedia', ztmp(:,:,1:16) ) zbg(1:16) = glob_sum_vec( 'icedia', ztmp(:,:,1:16) )
...@@ -261,9 +266,9 @@ CONTAINS ...@@ -261,9 +266,9 @@ CONTAINS
frc_tembot = 0._wp frc_tembot = 0._wp
frc_sal = 0._wp frc_sal = 0._wp
! record initial ice volume, salt and temp ! record initial ice volume, salt and temp
vol_loc_ini(:,:) = rhoi * vt_i(:,:) + rhos * vt_s(:,:) ! ice/snow volume (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) tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:) ! ice/snow heat content (J)
sal_loc_ini(:,:) = rhoi * st_i(:,:) ! ice salt content (pss*kg/m2) sal_loc_ini(:,:) = rhoi * st_i(:,:) ! ice salt content (pss*kg/m2)
ENDIF ENDIF
! !
ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
......
src/ICE/icedyn.F90
View file @ 3efbcb17
...@@ -93,8 +93,8 @@ CONTAINS ...@@ -93,8 +93,8 @@ CONTAINS
! !
! retrieve thickness from volume for landfast param. and UMx advection scheme ! retrieve thickness from volume for landfast param. and UMx advection scheme
WHERE( a_i(:,:,:) >= epsi20 ) WHERE( a_i(:,:,:) >= epsi20 )
h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:) h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:)
h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:) h_s(:,:,:) = v_s(:,:,:) / a_i(:,:,:)
ELSEWHERE ELSEWHERE
h_i(:,:,:) = 0._wp h_i(:,:,:) = 0._wp
h_s(:,:,:) = 0._wp h_s(:,:,:) = 0._wp
...@@ -162,15 +162,12 @@ CONTAINS ...@@ -162,15 +162,12 @@ CONTAINS
CASE ( np_dynADV1D , np_dynADV2D ) CASE ( np_dynADV1D , np_dynADV2D )
ALLOCATE( zdivu_i(jpi,jpj) ) ALLOCATE( zdivu_i(A2D(0)) )
DO_2D( 0, 0, 0, 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) & 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) & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
END_2D END_2D
CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1.0_wp )
! output
CALL iom_put( 'icediv' , zdivu_i ) CALL iom_put( 'icediv' , zdivu_i )
DEALLOCATE( zdivu_i ) DEALLOCATE( zdivu_i )
END SELECT END SELECT
......
src/ICE/icedyn_adv.F90
View file @ 3efbcb17
...@@ -41,6 +41,8 @@ MODULE icedyn_adv ...@@ -41,6 +41,8 @@ MODULE icedyn_adv
! ** namelist (namdyn_adv) ** ! ** namelist (namdyn_adv) **
INTEGER :: nn_UMx ! order of the UMx advection scheme INTEGER :: nn_UMx ! order of the UMx advection scheme
! !
!! * Substitutions
# include "do_loop_substitute.h90"
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018) !! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: icedyn_adv.F90 13472 2020-09-16 13:05:19Z smasson $ !! $Id: icedyn_adv.F90 13472 2020-09-16 13:05:19Z smasson $
...@@ -92,11 +94,11 @@ CONTAINS ...@@ -92,11 +94,11 @@ CONTAINS
!------------ !------------
! diagnostics ! 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_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 (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), 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(:,:,:) - sv_i_b(:,:,:) , 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 (:,:,:) - v_i_b (:,:,:) , 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 (:,:,:) - v_s_b (:,:,:) , 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('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('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) 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 ...@@ -89,7 +89,7 @@ CONTAINS
INTEGER :: icycle ! number of sub-timestep for the advection INTEGER :: icycle ! number of sub-timestep for the advection
REAL(wp) :: zdt, z1_dt ! - - REAL(wp) :: zdt, z1_dt ! - -
REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 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) :: zudy, zvdx
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max, zs_i, zsi_max 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 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: ze_i, zei_max
...@@ -100,7 +100,7 @@ CONTAINS ...@@ -100,7 +100,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es
REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei
!! diagnostics !! 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' IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
...@@ -129,8 +129,7 @@ CONTAINS ...@@ -129,8 +129,7 @@ CONTAINS
END DO END DO
CALL icemax4D( ze_i , zei_max ) CALL icemax4D( ze_i , zei_max )
CALL icemax4D( ze_s , zes_max ) CALL icemax4D( ze_s , zes_max )
CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp, zes_max, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp )
! !
! !
! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
...@@ -155,14 +154,14 @@ CONTAINS ...@@ -155,14 +154,14 @@ CONTAINS
DO jt = 1, icycle DO jt = 1, icycle
! diagnostics ! diagnostics
zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & zdiag_adv_mass(:,:) = SUM( pv_i (A2D(0),:) , dim=3 ) * rhoi + SUM( pv_s (A2D(0),:) , dim=3 ) * rhos &
& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow
zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi zdiag_adv_salt(:,:) = SUM( psv_i(A2D(0),:) , dim=3 ) * rhoi
zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(A2D(0),1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , 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) ! record at_i before advection (for open water)
zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) zati1(:,:) = SUM( pa_i(A2D(0),:), dim=3 )
! --- transported fields --- ! ! --- transported fields --- !
DO jl = 1, jpl DO jl = 1, jpl
...@@ -275,8 +274,8 @@ CONTAINS ...@@ -275,8 +274,8 @@ CONTAINS
& , z0oi , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp & ! ice age & , 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 ) & , 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 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 ) & , 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 & , 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 ) & , sxxe , 'T', 1._wp, syye , 'T', 1._wp, sxye , 'T', 1._wp )
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
IF( ln_pnd_lids ) THEN IF( ln_pnd_lids ) THEN
...@@ -317,7 +316,7 @@ CONTAINS ...@@ -317,7 +316,7 @@ CONTAINS
END DO END DO
! !
! derive open water from ice concentration ! 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 ) DO_2D( 0, 0, 0, 0 )
pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water 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 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
...@@ -325,13 +324,13 @@ CONTAINS ...@@ -325,13 +324,13 @@ CONTAINS
CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1.0_wp ) CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1.0_wp )
! !
! --- diagnostics --- ! ! --- diagnostics --- !
diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 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(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & & + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow &
& - zdiag_adv_mass(:,:) ) * z1_dt & - 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 & - zdiag_adv_salt(:,:) ) * z1_dt
diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , 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(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & & - SUM(SUM( pe_s(A2D(0),1:nlay_s,:) , dim=4 ), dim=3 ) &
& - zdiag_adv_heat(:,:) ) * z1_dt & - zdiag_adv_heat(:,:) ) * z1_dt
! !
! --- Ensure non-negative fields --- ! ! --- Ensure non-negative fields --- !
......
src/ICE/icedyn_adv_umx.F90
View file @ 3efbcb17
...@@ -104,7 +104,7 @@ CONTAINS ...@@ -104,7 +104,7 @@ CONTAINS
! !
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs
!! diagnostics !! 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' IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
...@@ -133,8 +133,7 @@ CONTAINS ...@@ -133,8 +133,7 @@ CONTAINS
END DO END DO
CALL icemax4D( ze_i , zei_max ) CALL icemax4D( ze_i , zei_max )
CALL icemax4D( ze_s , zes_max ) CALL icemax4D( ze_s , zes_max )
CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp ) CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp, zes_max, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp )
! !
! !
! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
...@@ -182,11 +181,11 @@ CONTAINS ...@@ -182,11 +181,11 @@ CONTAINS
DO jt = 1, icycle DO jt = 1, icycle
! diagnostics ! diagnostics
zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & zdiag_adv_mass(:,:) = SUM( pv_i (A2D(0),:) , dim=3 ) * rhoi + SUM( pv_s (A2D(0),:) , dim=3 ) * rhos &
& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow
zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi zdiag_adv_salt(:,:) = SUM( psv_i(A2D(0),:) , dim=3 ) * rhoi
zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(A2D(0),1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , 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) ! record at_i before advection (for open water)
zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
...@@ -370,8 +369,7 @@ CONTAINS ...@@ -370,8 +369,7 @@ CONTAINS
ELSE 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 ) 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 ENDIF
CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp ) CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp, pe_s, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp )
! !
!== Open water area ==! !== Open water area ==!
zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
...@@ -382,13 +380,13 @@ CONTAINS ...@@ -382,13 +380,13 @@ CONTAINS
CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp ) CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp )
! !
! --- diagnostics --- ! ! --- diagnostics --- !
diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 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(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & & + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow &
& - zdiag_adv_mass(:,:) ) * z1_dt & - 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 & - zdiag_adv_salt(:,:) ) * z1_dt
diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , 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(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & & - SUM(SUM( pe_s(A2D(0),1:nlay_s,:) , dim=4 ), dim=3 ) &
& - zdiag_adv_heat(:,:) ) * z1_dt & - zdiag_adv_heat(:,:) ) * z1_dt
! !
! --- Ensure non-negative fields and in-bound thicknesses --- ! ! --- Ensure non-negative fields and in-bound thicknesses --- !
......
src/ICE/icedyn_rdgrft.F90
View file @ 3efbcb17
...@@ -174,7 +174,7 @@ CONTAINS ...@@ -174,7 +174,7 @@ CONTAINS
! !
npti = 0 ; nptidx(:) = 0 npti = 0 ; nptidx(:) = 0
ipti = 0 ; iptidx(:) = 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 IF ( at_i(ji,jj) > epsi10 ) THEN
npti = npti + 1 npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji nptidx( npti ) = (jj - 1) * jpi + ji
...@@ -280,8 +280,15 @@ CONTAINS ...@@ -280,8 +280,15 @@ CONTAINS
CALL ice_dyn_1d2d( 2 ) ! --- Move to 2D arrays --- ! CALL ice_dyn_1d2d( 2 ) ! --- Move to 2D arrays --- !
ENDIF 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 ! controls
IF( sn_cfctl%l_prtctl ) CALL ice_prt3D('icedyn_rdgrft') ! prints IF( sn_cfctl%l_prtctl ) CALL ice_prt3D('icedyn_rdgrft') ! prints
...@@ -302,8 +309,9 @@ CONTAINS ...@@ -302,8 +309,9 @@ CONTAINS
!! ** Method : Compute the thickness distribution of the ice and open water !! ** Method : Compute the thickness distribution of the ice and open water
!! participating in ridging and of the resulting ridges. !! 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 INTEGER :: ji, jl ! dummy loop indices
REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar
...@@ -504,39 +512,43 @@ CONTAINS ...@@ -504,39 +512,43 @@ CONTAINS
END DO END DO
END DO END DO
! !
! 3) closing_gross IF( PRESENT( pclosing_net ) ) THEN
!----------------- !
! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. ! 3) closing_gross
! NOTE: 0 < aksum <= 1 !-----------------
WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.
ELSEWHERE ; closing_gross(1:npti) = 0._wp ! NOTE: 0 < aksum <= 1
END WHERE WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
ELSEWHERE ; closing_gross(1:npti) = 0._wp
! correction to closing rate if excessive ice removal END WHERE
!----------------------------------------------------
! Reduce the closing rate if more than 100% of any ice category would be removed ! correction to closing rate if excessive ice removal
! Reduce the opening rate in proportion !----------------------------------------------------
DO jl = 1, jpl ! 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 DO ji = 1, npti
zfac = apartf(ji,jl) * closing_gross(ji) * rDt_ice zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice
IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN IF( zfac < 0._wp ) THEN ! would lead to negative ato_i
closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice 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 ENDIF
END DO END DO
END DO !
ENDIF
! 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
! !
END SUBROUTINE rdgrft_prep END SUBROUTINE rdgrft_prep
...@@ -861,7 +873,8 @@ CONTAINS ...@@ -861,7 +873,8 @@ CONTAINS
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER :: ji, jj, jl ! dummy loop indices INTEGER :: ji, jj, jl ! dummy loop indices
REAL(wp) :: z1_3 ! local scalars 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 LOGICAL :: ln_str_R75
REAL(wp) :: zhi, zcp REAL(wp) :: zhi, zcp
...@@ -886,7 +899,7 @@ CONTAINS ...@@ -886,7 +899,7 @@ CONTAINS
! !
! Identify grid cells with ice ! Identify grid cells with ice
npti = 0 ; nptidx(:) = 0 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 IF ( at_i(ji,jj) > epsi10 ) THEN
npti = npti + 1 npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji nptidx( npti ) = (jj - 1) * jpi + ji
...@@ -899,7 +912,7 @@ CONTAINS ...@@ -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), ato_i_1d(1:npti) , ato_i )
CALL tab_2d_1d( npti, nptidx(1:npti), zstrength(1:npti) , strength ) 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 zaksum(1:npti) = apartf(1:npti,0) !clem: aksum should be defined in the header => local to module
DO jl = 1, jpl DO jl = 1, jpl
...@@ -956,12 +969,20 @@ CONTAINS ...@@ -956,12 +969,20 @@ CONTAINS
CALL tab_1d_2d( npti, nptidx(1:npti), zstrength(1:npti), strength ) CALL tab_1d_2d( npti, nptidx(1:npti), zstrength(1:npti), strength )
! !
ENDIF 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 ==! 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 ==! 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 END SELECT
! !
......
src/ICE/icedyn_rhg_eap.F90
View file @ 3efbcb17
This diff is collapsed.
src/ICE/icedyn_rhg_evp.F90
View file @ 3efbcb17
This diff is collapsed.
src/ICE/icedyn_rhg_vp.F90
View file @ 3efbcb17
...@@ -332,7 +332,10 @@ CONTAINS ...@@ -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) 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 ! 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)) ! 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) zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
...@@ -369,7 +372,10 @@ CONTAINS ...@@ -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) 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 ! 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)) ! 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) 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 ...@@ -37,6 +37,7 @@ MODULE iceistate
USE lib_mpp ! MPP library USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero) USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE fldread ! read input fields USE fldread ! read input fields
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
# if defined key_agrif # if defined key_agrif
USE agrif_oce USE agrif_oce
...@@ -113,11 +114,11 @@ CONTAINS ...@@ -113,11 +114,11 @@ CONTAINS
INTEGER :: ji, jj, jk, jl ! dummy loop indices INTEGER :: ji, jj, jk, jl ! dummy loop indices
REAL(wp) :: ztmelts, zsshadj, area REAL(wp) :: ztmelts, zsshadj, area
INTEGER , DIMENSION(4) :: itest INTEGER , DIMENSION(4) :: itest
REAL(wp), DIMENSION(jpi,jpj) :: z2d
REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator 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(A2D(0)) :: zmsk ! ice indicator
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(A2D(0)) :: zht_i_ini, zat_i_ini, ztm_s_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)) :: 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(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 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 ...@@ -135,53 +136,59 @@ CONTAINS
CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)
! !
! surface temperature and conductivity
DO jl = 1, jpl DO jl = 1, jpl
t_su (:,:,jl) = rt0 * tmask(:,:,1) ! temp at the surface ! == reduced arrays == !
cnd_ice(:,:,jl) = 0._wp ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T) DO_2D( 0, 0, 0, 0 )
END DO !
! cnd_ice(ji,jj,jl) = 0._wp ! conductivity at the ice top
! ice and snw temperatures !
DO jl = 1, jpl tn_ice(ji,jj,jl) = t_i (ji,jj,1,jl) ! temp for coupled runs
DO jk = 1, nlay_i t1_ice(ji,jj,jl) = t_i (ji,jj,1,jl) ! temp for coupled runs
t_i(:,:,jk,jl) = rt0 * tmask(:,:,1) !
END DO a_ip_eff(ji,jj,jl) = 0._wp ! melt pond effective fraction
DO jk = 1, nlay_s END_2D
t_s(:,:,jk,jl) = rt0 * tmask(:,:,1) !
END DO ! == full arrays == !
END DO DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! !
! specific temperatures for coupled runs ! heat contents
tn_ice (:,:,:) = t_i (:,:,1,:) DO jk = 1, nlay_i
t1_ice (:,:,:) = t_i (:,:,1,:) e_i(ji,jj,jk,jl) = 0._wp
t_i(ji,jj,jk,jl) = rt0 * tmask(ji,jj,1) ! ice temp
! heat contents ENDDO
e_i (:,:,:,:) = 0._wp DO jk = 1, nlay_s
e_s (:,:,:,:) = 0._wp e_s(ji,jj,jk,jl) = 0._wp
t_s(ji,jj,jk,jl) = rt0 * tmask(ji,jj,1) ! snw temp
! general fields ENDDO
a_i (:,:,:) = 0._wp !
v_i (:,:,:) = 0._wp ! general fields
v_s (:,:,:) = 0._wp a_i (ji,jj,jl) = 0._wp
sv_i(:,:,:) = 0._wp v_i (ji,jj,jl) = 0._wp
oa_i(:,:,:) = 0._wp v_s (ji,jj,jl) = 0._wp
! sv_i(ji,jj,jl) = 0._wp
h_i (:,:,:) = 0._wp oa_i(ji,jj,jl) = 0._wp
h_s (:,:,:) = 0._wp h_i (ji,jj,jl) = 0._wp
s_i (:,:,:) = 0._wp h_s (ji,jj,jl) = 0._wp
o_i (:,:,:) = 0._wp s_i (ji,jj,jl) = 0._wp
! o_i (ji,jj,jl) = 0._wp
! melt ponds t_su(ji,jj,jl) = rt0 * tmask(ji,jj,1)
a_ip (:,:,:) = 0._wp !
v_ip (:,:,:) = 0._wp ! melt ponds
v_il (:,:,:) = 0._wp a_ip(ji,jj,jl) = 0._wp
a_ip_eff (:,:,:) = 0._wp v_ip(ji,jj,jl) = 0._wp
h_ip (:,:,:) = 0._wp v_il(ji,jj,jl) = 0._wp
h_il (:,:,:) = 0._wp h_ip(ji,jj,jl) = 0._wp
h_il(ji,jj,jl) = 0._wp
!
END_2D
!
ENDDO
! !
! ice velocities ! ice velocities
u_ice (:,:) = 0._wp DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
v_ice (:,:) = 0._wp 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 ! 2) overwrite some of the fields with namelist parameters or netcdf file
...@@ -194,30 +201,30 @@ CONTAINS ...@@ -194,30 +201,30 @@ CONTAINS
! !---------------! ! !---------------!
IF( nn_iceini_file == 1 )THEN ! Read a file ! IF( nn_iceini_file == 1 )THEN ! Read a file !
! !---------------! ! !---------------!
WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp WHERE( ff_t(A2D(0)) >= 0._wp ) ; zmsk(:,:) = 1._wp
ELSEWHERE ; zswitch(:,:) = 0._wp ELSEWHERE ; zmsk(:,:) = 0._wp
END WHERE END WHERE
! !
CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
! !
! -- mandatory fields -- ! ! -- mandatory fields -- !
zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1) zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * smask0(:,:)
zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1) zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * smask0(:,:)
zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1) zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * smask0(:,:)
! -- optional fields -- ! ! -- optional fields -- !
! if fields do not exist then set them to the values present in the namelist (except for temperatures) ! if fields do not exist then set them to the values present in the namelist (except for temperatures)
! !
! ice salinity ! ice salinity
IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) & 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 ! temperatures
IF ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. & 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 & 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_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 * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 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 * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zmsk + rn_tms_ini_s * (1._wp - zmsk) ) * smask0(:,:)
ENDIF 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 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 ) & si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 )
...@@ -234,67 +241,60 @@ CONTAINS ...@@ -234,67 +241,60 @@ CONTAINS
! !
! pond concentration ! pond concentration
IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 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) & * si(jp_ati)%fnow(:,:,1)
! !
! pond depth ! pond depth
IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 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 ! pond lid depth
IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) & 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) zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * smask0(:,:)
ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1) ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * smask0(:,:)
zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1) zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * smask0(:,:)
ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1) ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * smask0(:,:)
zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1) zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * smask0(:,:)
zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1) zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * smask0(:,:)
zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1) * tmask(:,:,1) 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 ! ELSE ! Read namelist !
! !---------------! ! !---------------!
! no ice if (sst - Tfreez) >= thresold ! no ice if (sst - Tfreez) >= thresold
WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp WHERE( ( sst_m(A2D(0)) - (t_bo(A2D(0)) - rt0) ) * smask0(:,:) >= rn_thres_sst ) ; zmsk(:,:) = 0._wp
ELSEWHERE ; zswitch(:,:) = tmask(:,:,1) ELSEWHERE ; zmsk(:,:) = smask0(:,:)
END WHERE END WHERE
! !
! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
WHERE( ff_t(:,:) >= 0._wp ) WHERE( ff_t(A2D(0)) >= 0._wp )
zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:) zht_i_ini(:,:) = rn_hti_ini_n * zmsk(:,:)
zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) zht_s_ini(:,:) = rn_hts_ini_n * zmsk(:,:)
zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) zat_i_ini(:,:) = rn_ati_ini_n * zmsk(:,:)
zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) zsm_i_ini(:,:) = rn_smi_ini_n * zmsk(:,:)
ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) ztm_i_ini(:,:) = rn_tmi_ini_n * zmsk(:,:)
zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) zt_su_ini(:,:) = rn_tsu_ini_n * zmsk(:,:)
ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) ztm_s_ini(:,:) = rn_tms_ini_n * zmsk(:,:)
zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 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 * zswitch(:,:) zhpnd_ini(:,:) = rn_hpd_ini_n * zmsk(:,:)
zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:) zhlid_ini(:,:) = rn_hld_ini_n * zmsk(:,:)
ELSEWHERE ELSEWHERE
zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) zht_i_ini(:,:) = rn_hti_ini_s * zmsk(:,:)
zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) zht_s_ini(:,:) = rn_hts_ini_s * zmsk(:,:)
zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) zat_i_ini(:,:) = rn_ati_ini_s * zmsk(:,:)
zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) zsm_i_ini(:,:) = rn_smi_ini_s * zmsk(:,:)
ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) ztm_i_ini(:,:) = rn_tmi_ini_s * zmsk(:,:)
zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:) zt_su_ini(:,:) = rn_tsu_ini_s * zmsk(:,:)
ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:) ztm_s_ini(:,:) = rn_tms_ini_s * zmsk(:,:)
zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 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 * zswitch(:,:) zhpnd_ini(:,:) = rn_hpd_ini_s * zmsk(:,:)
zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:) zhlid_ini(:,:) = rn_hld_ini_s * zmsk(:,:)
END WHERE END WHERE
! !
ENDIF ENDIF
! make sure ponds = 0 if no ponds scheme ! make sure ponds = 0 if no ponds scheme
IF ( .NOT.ln_pnd ) THEN IF ( .NOT.ln_pnd ) THEN
zapnd_ini(:,:) = 0._wp zapnd_ini(:,:) = 0._wp
...@@ -311,7 +311,7 @@ CONTAINS ...@@ -311,7 +311,7 @@ CONTAINS
!----------------! !----------------!
! select ice covered grid points ! select ice covered grid points
npti = 0 ; nptidx(:) = 0 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 IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
npti = npti + 1 npti = npti + 1
nptidx(npti) = (jj - 1) * jpi + ji nptidx(npti) = (jj - 1) * jpi + ji
...@@ -363,6 +363,16 @@ CONTAINS ...@@ -363,6 +363,16 @@ CONTAINS
DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
& zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_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 ! calculate extensive and intensive variables
CALL ice_var_salprof ! for sz_i CALL ice_var_salprof ! for sz_i
DO jl = 1, jpl DO jl = 1, jpl
...@@ -398,15 +408,17 @@ CONTAINS ...@@ -398,15 +408,17 @@ CONTAINS
ENDIF ENDIF
#endif #endif
! Melt ponds ! Melt ponds
WHERE( a_i > epsi10 ) ; a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) WHERE( a_i(A2D(0),:) > epsi10 ) ; a_ip_eff(:,:,:) = a_ip(A2D(0),:) / a_i(A2D(0),:)
ELSEWHERE ; a_ip_eff(:,:,:) = 0._wp ELSEWHERE ; a_ip_eff(:,:,:) = 0._wp
END WHERE END WHERE
v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:) v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:)
! specific temperatures for coupled runs ! specific temperatures for coupled runs
tn_ice(:,:,:) = t_su(:,:,:) DO_2D( 0, 0, 0, 0 )
t1_ice(:,:,:) = t_i (:,:,1,:) 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 ! ice concentration should not exceed amax
at_i(:,:) = SUM( a_i, dim=3 ) at_i(:,:) = SUM( a_i, dim=3 )
...@@ -546,8 +558,8 @@ CONTAINS ...@@ -546,8 +558,8 @@ CONTAINS
ENDIF ENDIF
! !
DO ifpr = 1, jpfldi DO ifpr = 1, jpfldi
ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) ) ALLOCATE( si(ifpr)%fnow(A2D(0),1) )
IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(A2D(0),1,2) )
END DO END DO
! !
! fill si with slf_i and control print ! fill si with slf_i and control print
......
src/ICE/iceitd.F90
View file @ 3efbcb17
...@@ -29,6 +29,7 @@ MODULE iceitd ...@@ -29,6 +29,7 @@ MODULE iceitd
USE lib_fortran ! fortran utilities (glob_sum + no signed zero) USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE prtctl ! Print control USE prtctl ! Print control
USE timing ! Timing USE timing ! Timing
USE lbclnk ! lateral boundary conditions (or mpp links)
IMPLICIT NONE IMPLICIT NONE
PRIVATE PRIVATE
...@@ -100,7 +101,7 @@ CONTAINS ...@@ -100,7 +101,7 @@ CONTAINS
at_i(:,:) = SUM( a_i, dim=3 ) at_i(:,:) = SUM( a_i, dim=3 )
! !
npti = 0 ; nptidx(:) = 0 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 IF ( at_i(ji,jj) > epsi10 ) THEN
npti = npti + 1 npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji nptidx( npti ) = (jj - 1) * jpi + ji
...@@ -327,6 +328,9 @@ CONTAINS ...@@ -327,6 +328,9 @@ CONTAINS
! !
ENDIF 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_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_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') IF( ln_timing ) CALL timing_stop ('iceitd_rem')
...@@ -626,7 +630,7 @@ CONTAINS ...@@ -626,7 +630,7 @@ CONTAINS
DO jl = 1, jpl-1 ! identify thicknesses that are too big DO jl = 1, jpl-1 ! identify thicknesses that are too big
! !--------------------------------------- ! !---------------------------------------
npti = 0 ; nptidx(:) = 0 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 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 npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji nptidx( npti ) = (jj - 1) * jpi + ji
...@@ -662,7 +666,7 @@ CONTAINS ...@@ -662,7 +666,7 @@ CONTAINS
DO jl = jpl-1, 1, -1 ! Identify thicknesses that are too small DO jl = jpl-1, 1, -1 ! Identify thicknesses that are too small
! !----------------------------------------- ! !-----------------------------------------
npti = 0 ; nptidx(:) = 0 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 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 npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji nptidx( npti ) = (jj - 1) * jpi + ji
...@@ -687,6 +691,14 @@ CONTAINS ...@@ -687,6 +691,14 @@ CONTAINS
! !
END DO 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_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_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') IF( ln_timing ) CALL timing_stop ('iceitd_reb')
......
src/ICE/icesbc.F90
View file @ 3efbcb17
...@@ -58,7 +58,7 @@ CONTAINS ...@@ -58,7 +58,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: utau_ice, vtau_ice ! air-ice stress [N/m2] REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: utau_ice, vtau_ice ! air-ice stress [N/m2]
!! !!
INTEGER :: ji, jj ! dummy loop index 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') IF( ln_timing ) CALL timing_start('icesbc')
...@@ -70,25 +70,38 @@ CONTAINS ...@@ -70,25 +70,38 @@ CONTAINS
ENDIF ENDIF
! !
SELECT CASE( ksbc ) SELECT CASE( ksbc )
CASE( jp_usr ) ; CALL usrdef_sbc_ice_tau( kt ) ! user defined formulation !
CASE( jp_blk ) CASE( jp_usr ) !--- User defined formulation
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... CALL usrdef_sbc_ice_tau( kt )
& sf(jp_slp )%fnow(:,:,1), u_ice, v_ice, tm_su , & ! inputs !
& putaui = utau_ice, pvtaui = vtau_ice ) ! outputs CASE( jp_blk ) !--- Forced formulation
! 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 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 END SELECT
! !
IF( ln_mixcpl) THEN ! Case of a mixed Bulk/Coupled formulation IF( ln_mixcpl) THEN !--- Case of a mixed Bulk/Coupled formulation
CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice ) !
CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
!
DO_2D( 0, 0, 0, 0 ) 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) ) 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) ) vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
END_2D END_2D
CALL lbc_lnk( 'icesbc', utau_ice, 'U', -1.0_wp, vtau_ice, 'V', -1.0_wp ) !
ENDIF ENDIF
! !
CALL lbc_lnk( 'icesbc', utau_ice, 'U', -1.0_wp, vtau_ice, 'V', -1.0_wp )
!
IF( ln_timing ) CALL timing_stop('icesbc') IF( ln_timing ) CALL timing_stop('icesbc')
! !
END SUBROUTINE ice_sbc_tau END SUBROUTINE ice_sbc_tau
...@@ -129,25 +142,49 @@ CONTAINS ...@@ -129,25 +142,49 @@ CONTAINS
WRITE(numout,*)'~~~~~~~~~~~~~~~' WRITE(numout,*)'~~~~~~~~~~~~~~~'
ENDIF ENDIF
! !== ice albedo ==! ! !== 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 ==! SELECT CASE( ksbc ) !== fluxes over sea ice ==!
! !
CASE( jp_usr ) !--- user defined formulation CASE( jp_usr ) !--- user defined formulation
!
CALL usrdef_sbc_ice_flx( kt, h_s, h_i ) CALL usrdef_sbc_ice_flx( kt, h_s, h_i )
!
CASE( jp_blk, jp_abl ) !--- bulk formulation & ABL formulation 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... CALL blk_ice_2( t_su(A2D(0),:), h_s(A2D(0),:), h_i(A2D(0),:), & ! <<== in
& sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), & & alb_ice(:,:,:), theta_air_zt(:,:), q_air_zt(:,:), & ! <<== in
& sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) ) & sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), & ! <<== in
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 ) & sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) ) ! <<== in
IF( nn_flxdist /= -1 ) CALL ice_flx_dist ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist ) !
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) ! ! compute conduction flux and surface temperature (as in Jules surface module)
IF( ln_cndflx .AND. .NOT.ln_cndemulate ) & IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN
& CALL blk_ice_qcn ( ln_virtual_itd, t_su, t_bo, h_s, h_i ) 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 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 END SELECT
! !== some fluxes at the ice-ocean interface and in the leads ! !== some fluxes at the ice-ocean interface and in the leads
CALL ice_flx_other CALL ice_flx_other
...@@ -157,7 +194,8 @@ CONTAINS ...@@ -157,7 +194,8 @@ CONTAINS
END SUBROUTINE ice_sbc_flx 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 *** !! *** ROUTINE ice_flx_dist ***
!! !!
...@@ -173,18 +211,20 @@ CONTAINS ...@@ -173,18 +211,20 @@ CONTAINS
!! using T-ice and albedo sensitivity !! using T-ice and albedo sensitivity
!! = 2 Redistribute a single flux over categories !! = 2 Redistribute a single flux over categories
!!------------------------------------------------------------------- !!-------------------------------------------------------------------
INTEGER , INTENT(in ) :: k_flxdist ! redistributor INTEGER , INTENT(in ) :: k_flxdist ! redistributor
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptn_ice ! ice surface temperature REAL(wp), DIMENSION(A2D(0)) , INTENT(in ) :: pat_i ! ice concentration
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb_ice ! ice albedo REAL(wp), DIMENSION(A2D(0),jpl), INTENT(in ) :: pa_i ! ice concentration
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqns_ice ! non solar flux REAL(wp), DIMENSION(A2D(0),jpl), INTENT(in ) :: ptn_ice ! ice surface temperature
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqsr_ice ! net solar flux REAL(wp), DIMENSION(A2D(0),jpl), INTENT(in ) :: palb_ice ! ice albedo
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdqn_ice ! non solar flux sensitivity REAL(wp), DIMENSION(A2D(0),jpl), INTENT(inout) :: pqns_ice ! non solar flux
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pevap_ice ! sublimation REAL(wp), DIMENSION(A2D(0),jpl), INTENT(inout) :: pqsr_ice ! net solar flux
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdevap_ice ! sublimation sensitivity 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 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_qsr_m ! Mean solar heat flux over all categories
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_qns_m ! Mean non 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 ...@@ -195,7 +235,7 @@ CONTAINS
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztem_m ! Mean temperature over all categories 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 ELSEWHERE ; z1_at_i(:,:) = 0._wp
END WHERE END WHERE
...@@ -203,13 +243,13 @@ CONTAINS ...@@ -203,13 +243,13 @@ CONTAINS
! !
CASE( 0 , 1 ) 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_qns_m (:,:) = SUM( pa_i(:,:,:) * pqns_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_qsr_m (:,:) = SUM( a_i(:,:,:) * pqsr_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) z_qsr_m (:,:) = SUM( pa_i(:,:,:) * pqsr_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_dqn_m (:,:) = SUM( a_i(:,:,:) * pdqn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) z_dqn_m (:,:) = SUM( pa_i(:,:,:) * pdqn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) z_evap_m (:,:) = SUM( pa_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:) z_devap_m(:,:) = SUM( pa_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
DO jl = 1, jpl DO jl = 1, jpl
pqns_ice (:,:,jl) = z_qns_m (:,:) pqns_ice (:,:,jl) = z_qns_m (:,:)
pqsr_ice (:,:,jl) = z_qsr_m (:,:) pqsr_ice (:,:,jl) = z_qsr_m (:,:)
...@@ -226,10 +266,10 @@ CONTAINS ...@@ -226,10 +266,10 @@ CONTAINS
! !
CASE( 1 , 2 ) 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(:,:) zalb_m(:,:) = SUM( pa_i(:,:,:) * palb_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
ztem_m(:,:) = SUM( a_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) ztem_m(:,:) = SUM( pa_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
DO jl = 1, jpl DO jl = 1, jpl
pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) ) 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(:,:) ) pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
...@@ -257,7 +297,7 @@ CONTAINS ...@@ -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) :: 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 :: 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), 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 ! computation of friction velocity at T points
...@@ -268,36 +308,33 @@ CONTAINS ...@@ -268,36 +308,33 @@ CONTAINS
zv_io = v_ice(ji ,jj ) - ssv_m(ji ,jj ) zv_io = v_ice(ji ,jj ) - ssv_m(ji ,jj )
zv_iom1 = v_ice(ji ,jj-1) - ssv_m(ji ,jj-1) 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) ) + & 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) ) ) & ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )
END_2D END_2D
ELSE ! if no ice dynamics => transfer directly the atmospheric stress to the ocean ELSE ! if no ice dynamics => transfer directly the atmospheric stress to the ocean
DO_2D( 0, 0, 0, 0 ) DO_2D( 0, 0, 0, 0 )
zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp * & zfric(ji,jj) = r1_rho0 * SQRT( utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj) ) * smask0(ji,jj)
& ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & zvel (ji,jj) = 0._wp
& + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
zvel(ji,jj) = 0._wp
END_2D END_2D
ENDIF 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 ! Partial computation of forcing for the thermodynamic sea ice model
!--------------------------------------------------------------------! !--------------------------------------------------------------------!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! needed for qlead DO_2D( 0, 0, 0, 0 ) ! needed for qlead
rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 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) --- ! ! --- 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) ) * qsr_oce(ji,jj) * frq_m(ji,jj) + &
& ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(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) --- ! ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- !
! (mostly<0 but >0 if supercooling) ! (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 = 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_neg = MIN( zqfr , 0._wp ) ! only < 0
zqfr_pos = MAX( zqfr , 0._wp ) ! only > 0 zqfr_pos = MAX( zqfr , 0._wp ) ! only > 0
! --- Sensible ocean-to-ice heat flux (W/m2) --- ! ! --- Sensible ocean-to-ice heat flux (W/m2) --- !
! (mostly>0 but <0 if supercooling) ! (mostly>0 but <0 if supercooling)
......
src/ICE/icestp.F90
View file @ 3efbcb17
...@@ -77,6 +77,7 @@ MODULE icestp ...@@ -77,6 +77,7 @@ MODULE icestp
USE lib_fortran ! fortran utilities (glob_sum + no signed zero) USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE timing ! Timing USE timing ! Timing
USE prtctl ! Print control USE prtctl ! Print control
USE lbclnk ! lateral boundary conditions (or mpp links)
IMPLICIT NONE IMPLICIT NONE
PRIVATE PRIVATE
...@@ -109,7 +110,7 @@ CONTAINS ...@@ -109,7 +110,7 @@ CONTAINS
!! - save the outputs !! - save the outputs
!! - save the outputs for restart when necessary !! - 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: !! - update all sbc variables below sea-ice:
!! utau, vtau, taum, wndm, qns , qsr, emp , sfx !! utau, vtau, taum, wndm, qns , qsr, emp , sfx
!!--------------------------------------------------------------------- !!---------------------------------------------------------------------
...@@ -117,7 +118,7 @@ CONTAINS ...@@ -117,7 +118,7 @@ CONTAINS
INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices
INTEGER, INTENT(in) :: ksbc ! flux formulation (user defined, bulk, or Pure Coupled) 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') IF( ln_timing ) CALL timing_start('icestp')
...@@ -128,9 +129,12 @@ CONTAINS ...@@ -128,9 +129,12 @@ CONTAINS
! !
kt_ice = kt ! -- Ice model time step kt_ice = kt ! -- Ice model time step
! !
u_oce(:,:) = ssu_m(:,:) ! -- mean surface ocean current DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! -- mean surface ocean current
v_oce(:,:) = ssv_m(:,:) 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) 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) ) t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
! !
...@@ -152,6 +156,7 @@ CONTAINS ...@@ -152,6 +156,7 @@ CONTAINS
! utau_ice, vtau_ice = surface ice stress [N/m2] ! utau_ice, vtau_ice = surface ice stress [N/m2]
!------------------------------------------------! !------------------------------------------------!
CALL ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice ) CALL ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice )
! => clem: here utau_ice and vtau_ice are defined everywhere
!-------------------------------------! !-------------------------------------!
! --- ice dynamics and advection --- ! ! --- ice dynamics and advection --- !
!-------------------------------------! !-------------------------------------!
...@@ -161,6 +166,11 @@ CONTAINS ...@@ -161,6 +166,11 @@ CONTAINS
IF( ln_icedyn .AND. .NOT.ln_c1d ) & IF( ln_icedyn .AND. .NOT.ln_c1d ) &
& CALL ice_dyn( kt, Kmm ) ! -- Ice dynamics & 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 CALL diag_trends( 1 ) ! record dyn trends
! !
! !== lateral boundary conditions ==! ! !== lateral boundary conditions ==!
...@@ -170,7 +180,8 @@ CONTAINS ...@@ -170,7 +180,8 @@ CONTAINS
CALL ice_var_glo2eqv ! h_i and h_s for ice albedo calculation CALL ice_var_glo2eqv ! h_i and h_s for ice albedo calculation
CALL ice_var_agg(1) ! at_i for coupling CALL ice_var_agg(1) ! at_i for coupling
CALL store_fields ! Store now ice values 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 --- ! ! --- Thermodynamical coupling with the atmosphere --- !
!------------------------------------------------------! !------------------------------------------------------!
...@@ -184,11 +195,15 @@ CONTAINS ...@@ -184,11 +195,15 @@ CONTAINS
! qemp_oce, qemp_ice, = sensible heat (associated with evap & precip) [W/m2] ! qemp_oce, qemp_ice, = sensible heat (associated with evap & precip) [W/m2]
! qprec_ice, qevap_ice ! 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 ) CALL ice_sbc_flx( kt, ksbc )
!----------------------------! !----------------------------!
! --- ice thermodynamics --- ! ! --- ice thermodynamics --- !
!----------------------------! !----------------------------!
IF( ln_icethd ) CALL ice_thd( kt ) ! -- 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 diag_trends( 2 ) ! record thermo trends
CALL ice_var_glo2eqv ! necessary calls (at least for coupling) CALL ice_var_glo2eqv ! necessary calls (at least for coupling)
...@@ -288,7 +303,7 @@ CONTAINS ...@@ -288,7 +303,7 @@ CONTAINS
CALL ice_drift_init ! initialization for diags of conservation CALL ice_drift_init ! initialization for diags of conservation
! !
fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction 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 IF( ln_rstart ) THEN
CALL iom_close( numrir ) ! close input ice restart file CALL iom_close( numrir ) ! close input ice restart file
...@@ -369,26 +384,33 @@ CONTAINS ...@@ -369,26 +384,33 @@ CONTAINS
INTEGER :: ji, jj, jl ! dummy loop index INTEGER :: ji, jj, jl ! dummy loop index
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
! !
a_i_b (:,:,:) = a_i (:,:,:) ! ice area DO jl = 1, jpl
v_i_b (:,:,:) = v_i (:,:,:) ! ice volume DO_2D( 0, 0, 0, 0 )
v_s_b (:,:,:) = v_s (:,:,:) ! snow volume a_i_b (ji,jj,jl) = a_i (ji,jj,jl) ! ice area
v_ip_b(:,:,:) = v_ip(:,:,:) ! pond volume v_i_b (ji,jj,jl) = v_i (ji,jj,jl) ! ice volume
v_il_b(:,:,:) = v_il(:,:,:) ! pond lid volume v_s_b (ji,jj,jl) = v_s (ji,jj,jl) ! snow volume
sv_i_b(:,:,:) = sv_i(:,:,:) ! salt content v_ip_b(ji,jj,jl) = v_ip(ji,jj,jl) ! pond volume
e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy v_il_b(ji,jj,jl) = v_il(ji,jj,jl) ! pond lid volume
e_i_b (:,:,:,:) = e_i (:,:,:,:) ! ice thermal energy sv_i_b(ji,jj,jl) = sv_i(ji,jj,jl) ! salt content
WHERE( a_i_b(:,:,:) >= epsi20 ) IF( a_i_b(ji,jj,jl) >= epsi20 ) THEN
h_i_b(:,:,:) = v_i_b(:,:,:) / a_i_b(:,:,:) ! ice thickness h_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / a_i_b(ji,jj,jl) ! ice thickness
h_s_b(:,:,:) = v_s_b(:,:,:) / a_i_b(:,:,:) ! snw thickness h_s_b(ji,jj,jl) = v_s_b(ji,jj,jl) / a_i_b(ji,jj,jl) ! snw thickness
ELSEWHERE ELSE
h_i_b(:,:,:) = 0._wp h_i_b(ji,jj,jl) = 0._wp
h_s_b(:,:,:) = 0._wp h_s_b(ji,jj,jl) = 0._wp
END WHERE ENDIF
! e_s_b (ji,jj,:,jl) = e_s (ji,jj,:,jl) ! snow thermal energy
! ice velocities & total concentration 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 ) 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 END SUBROUTINE store_fields
...@@ -402,20 +424,23 @@ CONTAINS ...@@ -402,20 +424,23 @@ CONTAINS
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER :: ji, jj, jl ! dummy loop index 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 (ji,jj) = 0._wp ;
sfx_bri(ji,jj) = 0._wp ; sfx_lam(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_sni(ji,jj) = 0._wp ; sfx_opw(ji,jj) = 0._wp
sfx_bog(ji,jj) = 0._wp ; sfx_dyn(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_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_snw(ji,jj) = 0._wp ; wfx_ice(ji,jj) = 0._wp
wfx_sni(ji,jj) = 0._wp ; wfx_opw(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_bog(ji,jj) = 0._wp ; wfx_dyn(ji,jj) = 0._wp
wfx_bom(ji,jj) = 0._wp ; wfx_sum(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_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_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 wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp
...@@ -426,7 +451,7 @@ CONTAINS ...@@ -426,7 +451,7 @@ CONTAINS
hfx_snw(ji,jj) = 0._wp ; hfx_opw(ji,jj) = 0._wp 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_bog(ji,jj) = 0._wp ; hfx_dyn(ji,jj) = 0._wp
hfx_bom(ji,jj) = 0._wp ; hfx_sum(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_spr(ji,jj) = 0._wp ; hfx_dif(ji,jj) = 0._wp
hfx_err_dif(ji,jj) = 0._wp hfx_err_dif(ji,jj) = 0._wp
wfx_err_sub(ji,jj) = 0._wp wfx_err_sub(ji,jj) = 0._wp
...@@ -451,7 +476,7 @@ CONTAINS ...@@ -451,7 +476,7 @@ CONTAINS
END_2D END_2D
DO jl = 1, jpl DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) DO_2D( 0, 0, 0, 0 )
! SIMIP diagnostics ! SIMIP diagnostics
t_si (ji,jj,jl) = rt0 ! temp at the ice-snow interface t_si (ji,jj,jl) = rt0 ! temp at the ice-snow interface
qcn_ice_bot(ji,jj,jl) = 0._wp qcn_ice_bot(ji,jj,jl) = 0._wp
...@@ -477,22 +502,25 @@ CONTAINS ...@@ -477,22 +502,25 @@ CONTAINS
!! and outputs !! and outputs
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kn ! 1 = after dyn ; 2 = after thermo 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) ! --- trends of heat, salt, mass (used for conservation controls)
IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN
! !
diag_heat(:,:) = diag_heat(:,:) & DO_2D( 0, 0, 0, 0 )
& - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & diag_heat(ji,jj) = diag_heat(ji,jj) &
& - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice & - SUM(SUM( e_i (ji,jj,1:nlay_i,:) - e_i_b (ji,jj,1:nlay_i,:), dim=2 ) ) * r1_Dt_ice &
diag_sice(:,:) = diag_sice(:,:) & & - SUM(SUM( e_s (ji,jj,1:nlay_s,:) - e_s_b (ji,jj,1:nlay_s,:), dim=2 ) ) * r1_Dt_ice
& + SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice * rhoi diag_sice(ji,jj) = diag_sice(ji,jj) &
diag_vice(:,:) = diag_vice(:,:) & & + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * r1_Dt_ice * rhoi
& + SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhoi diag_vice(ji,jj) = diag_vice(ji,jj) &
diag_vsnw(:,:) = diag_vsnw(:,:) & & + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * r1_Dt_ice * rhoi
& + SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos diag_vsnw(ji,jj) = diag_vsnw(ji,jj) &
diag_vpnd(:,:) = diag_vpnd(:,:) & & + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * r1_Dt_ice * rhos
& + SUM( v_ip + v_il - v_ip_b - v_il_b , dim=3 ) * r1_Dt_ice * rhow 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 IF( kn == 2 ) CALL iom_put ( 'hfxdhc' , diag_heat ) ! output of heat trend
! !
...@@ -501,10 +529,12 @@ CONTAINS ...@@ -501,10 +529,12 @@ CONTAINS
! --- trends of concentration (used for simip outputs) ! --- trends of concentration (used for simip outputs)
IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 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 == 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 IF( kn == 2 ) CALL iom_put( 'afxtot' , diag_aice ) ! total trend
! !
ENDIF ENDIF
......
src/ICE/icetab.F90
View file @ 3efbcb17
...@@ -38,15 +38,24 @@ CONTAINS ...@@ -38,15 +38,24 @@ CONTAINS
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ndim1d ! 1d size INTEGER , INTENT(in ) :: ndim1d ! 1d size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 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 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 jl = 1, jpl
DO jn = 1, ndim1d DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
jjd = ( tab_ind(jn) - 1 ) / jpi + 1 jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
tab1d(jn,jl) = tab2d(jid,jjd,jl) tab1d(jn,jl) = tab2d(jid,jjd,jl)
END DO END DO
END DO END DO
...@@ -59,14 +68,23 @@ CONTAINS ...@@ -59,14 +68,23 @@ CONTAINS
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ndim1d ! 1d size INTEGER , INTENT(in ) :: ndim1d ! 1d size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 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 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 DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
jjd = ( tab_ind(jn) - 1 ) / jpi + 1 jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
tab1d( jn) = tab2d( jid, jjd) tab1d( jn) = tab2d( jid, jjd)
END DO END DO
END SUBROUTINE tab_2d_1d END SUBROUTINE tab_2d_1d
...@@ -79,14 +97,23 @@ CONTAINS ...@@ -79,14 +97,23 @@ CONTAINS
INTEGER , INTENT(in ) :: ndim1d ! 1D size INTEGER , INTENT(in ) :: ndim1d ! 1D size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in ) :: tab1d ! input 1D field 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 jl = 1, jpl
DO jn = 1, ndim1d DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
jjd = ( tab_ind(jn) - 1 ) / jpi + 1 jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
tab2d(jid,jjd,jl) = tab1d(jn,jl) tab2d(jid,jjd,jl) = tab1d(jn,jl)
END DO END DO
END DO END DO
...@@ -100,13 +127,22 @@ CONTAINS ...@@ -100,13 +127,22 @@ CONTAINS
INTEGER , INTENT(in ) :: ndim1d ! 1D size INTEGER , INTENT(in ) :: ndim1d ! 1D size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
REAL(wp), DIMENSION(ndim1d) , INTENT(in ) :: tab1d ! input 1D field 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 DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
jjd = ( tab_ind(jn) - 1 ) / jpi + 1 jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
tab2d(jid, jjd) = tab1d( jn) tab2d(jid, jjd) = tab1d( jn)
END DO END DO
END SUBROUTINE tab_1d_2d END SUBROUTINE tab_1d_2d
......
src/ICE/icethd.F90
View file @ 3efbcb17
...@@ -98,7 +98,7 @@ CONTAINS ...@@ -98,7 +98,7 @@ CONTAINS
! convergence tests ! convergence tests
IF( ln_zdf_chkcvg ) THEN 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 ztice_cvgerr = 0._wp ; ztice_cvgstp = 0._wp
ENDIF ENDIF
! !
...@@ -111,7 +111,7 @@ CONTAINS ...@@ -111,7 +111,7 @@ CONTAINS
! select ice covered grid points ! select ice covered grid points
npti = 0 ; nptidx(:) = 0 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 IF ( a_i(ji,jj,jl) > epsi10 ) THEN
npti = npti + 1 npti = npti + 1
nptidx(npti) = (jj - 1) * jpi + ji nptidx(npti) = (jj - 1) * jpi + ji
...@@ -158,6 +158,13 @@ CONTAINS ...@@ -158,6 +158,13 @@ CONTAINS
IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- ! 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 --- ! 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 --- ! CALL ice_cor( kt , 2 ) ! --- Corrections --- !
! !
......
src/ICE/icethd_do.F90
View file @ 3efbcb17
...@@ -28,7 +28,6 @@ MODULE icethd_do ...@@ -28,7 +28,6 @@ MODULE icethd_do
USE in_out_manager ! I/O manager USE in_out_manager ! I/O manager
USE lib_mpp ! MPP library USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero) USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE lbclnk ! lateral boundary conditions (or mpp links)
IMPLICIT NONE IMPLICIT NONE
PRIVATE PRIVATE
...@@ -105,7 +104,7 @@ CONTAINS ...@@ -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_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 ) 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 ! 1) Compute thickness, salinity, enthalpy, age, area and volume of new ice
!------------------------------------------------------------------------------! !------------------------------------------------------------------------------!
...@@ -113,7 +112,7 @@ CONTAINS ...@@ -113,7 +112,7 @@ CONTAINS
! Identify grid points where new ice forms ! Identify grid points where new ice forms
npti = 0 ; nptidx(:) = 0 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 IF ( qlead(ji,jj) < 0._wp ) THEN
npti = npti + 1 npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji nptidx( npti ) = (jj - 1) * jpi + ji
...@@ -325,6 +324,8 @@ CONTAINS ...@@ -325,6 +324,8 @@ CONTAINS
! !
ENDIF ! npti > 0 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_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) 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 ...@@ -372,8 +373,8 @@ CONTAINS
DO_2D( 0, 0, 0, 0 ) DO_2D( 0, 0, 0, 0 )
IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling
! -- Wind stress -- ! ! -- Wind stress -- !
ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) + utau_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp ztaux = utau_ice(ji,jj) * smask0(ji,jj)
ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + vtau_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp ztauy = vtau_ice(ji,jj) * smask0(ji,jj)
! Square root of wind stress ! Square root of wind stress
ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
...@@ -415,8 +416,6 @@ CONTAINS ...@@ -415,8 +416,6 @@ CONTAINS
! !
END_2D END_2D
! !
CALL lbc_lnk( 'icethd_frazil', fraz_frac, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp )
ENDIF ENDIF
END SUBROUTINE ice_thd_frazil END SUBROUTINE ice_thd_frazil
......