- Timestamp:
- 2017-09-27T12:09:10+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM
- Files:
-
- 14 edited
-
CONFIG/SHARED/namelist_ice_lim3_ref (modified) (1 diff)
-
NEMO/LIM_SRC_3/ice1D.F90 (modified) (2 diffs)
-
NEMO/LIM_SRC_3/icealb.F90 (modified) (3 diffs)
-
NEMO/LIM_SRC_3/icectl.F90 (modified) (5 diffs)
-
NEMO/LIM_SRC_3/icedyn_rdgrft.F90 (modified) (31 diffs)
-
NEMO/LIM_SRC_3/iceitd.F90 (modified) (26 diffs)
-
NEMO/LIM_SRC_3/icethd.F90 (modified) (5 diffs)
-
NEMO/LIM_SRC_3/icethd_da.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/icethd_dh.F90 (modified) (19 diffs)
-
NEMO/LIM_SRC_3/icethd_do.F90 (modified) (19 diffs)
-
NEMO/LIM_SRC_3/icethd_ent.F90 (modified) (4 diffs)
-
NEMO/LIM_SRC_3/icethd_sal.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/icethd_zdf.F90 (modified) (37 diffs)
-
NEMO/LIM_SRC_3/icevar.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r8531 r8565 68 68 rn_pstar = 2.0e+04 ! ice strength thickness parameter [N/m2] 69 69 rn_crhg = 20.0 ! ice strength conc. parameter (-) 70 ln_str_R75 = .false. ! ice strength param.: Rothrock_75 => P = Cf*coeff*integral(wr.h^2)71 rn_perdg = 17.0 ! ridging work divided by pot. energy change in ridging72 70 ! -- ice_rdgrft -- ! 73 71 rn_csrdg = 0.5 ! fraction of shearing energy contributing to ridging -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90
r8564 r8565 27 27 !: are the variables corresponding to 2d vectors 28 28 29 !! please, DOCTOR naming convention starting by i means LOCAL integer 30 !! ===>>> rename idxice as nidice or nidthd, or nidx_thd or midx ... ? 31 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: idxice !: selected points for ice thermo 32 INTEGER , PUBLIC :: nidx ! number of selected points 29 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: nptidx !: selected points for ice thermo 30 INTEGER , PUBLIC :: npti ! number of selected points 33 31 34 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qlead_1d … … 175 173 176 174 ii = 1 177 ALLOCATE( idxice(jpij) , &175 ALLOCATE( nptidx (jpij) , & 178 176 & qlead_1d (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d (jpij) , & 179 177 & fr1_i0_1d(jpij) , fr2_i0_1d (jpij) , qns_ice_1d(jpij) , & -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90
r8534 r8565 29 29 REAL(wp), PUBLIC, PARAMETER :: rn_alb_oce = 0.066 !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 30 30 31 REAL(wp) , PARAMETER :: rc1 = 0.05 ! snow thickness (only for nn_ice_alb=0)32 REAL(wp) , PARAMETER :: rc2 = 0.10 ! " "33 REAL(wp) , PARAMETER :: rcloud = 0.06 ! cloud effect on albedo (only-for nn_ice_alb=0)34 REAL(wp) , PARAMETER :: r1_c1 = 1. / rc135 REAL(wp) , PARAMETER :: r1_c2 = 1. / rc231 REAL(wp) , PARAMETER :: ppc1 = 0.05 ! snow thickness (only for nn_ice_alb=0) 32 REAL(wp) , PARAMETER :: ppc2 = 0.10 ! " " 33 REAL(wp) , PARAMETER :: ppcloud = 0.06 ! cloud effect on albedo (only-for nn_ice_alb=0) 34 REAL(wp) , PARAMETER :: pp1_c1 = 1. / ppc1 35 REAL(wp) , PARAMETER :: pp1_c2 = 1. / ppc2 36 36 ! 37 37 ! ** albedo namelist (namalb) … … 154 154 DO ji = 1, jpi 155 155 ! ! Freezing snow 156 ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > rc2157 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - rc1 ) ) )156 ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > ppc2 157 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - ppc1 ) ) ) 158 158 zalb_sf = ( 1._wp - zswitch ) * ( zalb_it(ji,jj,jl) & 159 & + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * r1_c1 ) &159 & + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * pp1_c1 ) & 160 160 & + zswitch * rn_alb_sdry 161 161 ! 162 162 ! ! Melting snow 163 ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > rc2164 zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - rc2 ) )165 zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * r1_c2 ) &163 ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > ppc2 164 zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - ppc2 ) ) 165 zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * pp1_c2 ) & 166 166 & + zswitch * rn_alb_smlt 167 167 ! … … 178 178 END DO 179 179 ! 180 pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud ! Oberhuber correction for overcast sky180 pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + ppcloud ! Oberhuber correction for overcast sky 181 181 ! 182 182 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
r8564 r8565 71 71 REAL(wp) :: zvtrp, zetrp 72 72 REAL(wp) :: zarea, zv_sill, zs_sill, zt_sill 73 REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt73 REAL(wp), PARAMETER :: ppconv = 1.e-9 ! convert W to GW and kg to Mt 74 74 !!------------------------------------------------------------------- 75 75 ! … … 79 79 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_ice_sub(:,:) + & 80 80 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:) & 81 & ) * e1e2t(:,:) ) * zconv81 & ) * e1e2t(:,:) ) * ppconv 82 82 ! 83 83 ! ! salt flux 84 84 pdiag_fs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 85 85 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) & 86 & ) * e1e2t(:,:) ) * zconv86 & ) * e1e2t(:,:) ) * ppconv 87 87 ! 88 88 ! ! heat flux 89 89 pdiag_ft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 90 90 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 91 & ) * e1e2t(:,:) ) * zconv92 93 pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv )94 95 pdiag_s = glob_sum( SUM( sv_i * rhoic , dim=3 ) * e1e2t * zconv )91 & ) * e1e2t(:,:) ) * ppconv 92 93 pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * ppconv ) 94 95 pdiag_s = glob_sum( SUM( sv_i * rhoic , dim=3 ) * e1e2t * ppconv ) 96 96 97 97 pdiag_t = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) & 98 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv98 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * ppconv 99 99 100 100 ELSEIF( icount == 1 ) THEN … … 104 104 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_ice_sub(:,:) + & 105 105 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:) & 106 & ) * e1e2t(:,:) ) * zconv - pdiag_fv106 & ) * e1e2t(:,:) ) * ppconv - pdiag_fv 107 107 108 108 ! salt flux 109 109 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 110 110 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) & 111 & ) * e1e2t(:,:) ) * zconv - pdiag_fs111 & ) * e1e2t(:,:) ) * ppconv - pdiag_fs 112 112 113 113 ! heat flux 114 114 zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 115 115 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 116 & ) * e1e2t(:,:) ) * zconv - pdiag_ft116 & ) * e1e2t(:,:) ) * ppconv - pdiag_ft 117 117 118 118 ! outputs 119 zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv &119 zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * ppconv & 120 120 & - pdiag_v ) * r1_rdtice - zfv ) * rday 121 121 122 zs = ( ( glob_sum( SUM( sv_i * rhoic , dim=3 ) * e1e2t ) * zconv &122 zs = ( ( glob_sum( SUM( sv_i * rhoic , dim=3 ) * e1e2t ) * ppconv & 123 123 & - pdiag_s ) * r1_rdtice + zfs ) * rday 124 124 125 125 zt = ( glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) & 126 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv &126 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * ppconv & 127 127 & - pdiag_t ) * r1_rdtice + zft 128 128 129 129 ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 130 zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t ) * zconv * rday131 zetrp = glob_sum( ( diag_trp_ei + diag_trp_es ) * e1e2t ) * zconv130 zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t ) * ppconv * rday 131 zetrp = glob_sum( ( diag_trp_ei + diag_trp_es ) * e1e2t ) * ppconv 132 132 133 133 zvmin = glob_min( v_i ) … … 136 136 137 137 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 138 zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2138 zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * ppconv ! in 1.e9 m2 139 139 zv_sill = zarea * 2.5e-5 140 140 zs_sill = zarea * 25.e-5 … … 176 176 REAL(wp) :: zhfx, zsfx, zvfx 177 177 REAL(wp) :: zarea, zv_sill, zs_sill, zt_sill 178 REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt178 REAL(wp), PARAMETER :: ppconv = 1.e-9 ! convert W to GW and kg to Mt 179 179 !!------------------------------------------------------------------- 180 180 181 181 ! water flux 182 zvfx = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday182 zvfx = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * ppconv * rday 183 183 184 184 ! salt flux 185 zsfx = glob_sum( ( sfx + diag_sice ) * e1e2t ) * zconv * rday185 zsfx = glob_sum( ( sfx + diag_sice ) * e1e2t ) * ppconv * rday 186 186 187 187 ! heat flux 188 188 zhfx = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es & 189 189 ! & - SUM( qevap_ice * a_i_b, dim=3 ) & !!clem: I think this line must be commented (but need check) 190 & ) * e1e2t ) * zconv190 & ) * e1e2t ) * ppconv 191 191 192 192 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 193 zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2193 zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * ppconv ! in 1.e9 m2 194 194 zv_sill = zarea * 2.5e-5 195 195 zs_sill = zarea * 25.e-5 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90
r8564 r8565 61 61 REAL(wp) :: rn_pstar ! determines ice strength, Hibler JPO79 62 62 REAL(wp) :: rn_crhg ! determines changes in ice strength 63 LOGICAL :: ln_str_R75 ! ice strength parameterization (Rothrock75)64 REAL(wp) :: rn_perdg ! ridging work divided by pot. energy change in ridging65 63 REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging 66 64 LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975)) … … 130 128 !! 131 129 INTEGER :: ji, jj, jk, jl ! dummy loop index 132 INTEGER :: niter, iterate_ridging! local integer133 INTEGER :: nidx2! local integer130 INTEGER :: iter, iterate_ridging ! local integer 131 INTEGER :: ipti ! local integer 134 132 REAL(wp) :: zfac ! local scalar 135 INTEGER , DIMENSION(jpij) :: i dxice2! compute ridge/raft or not133 INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not 136 134 REAL(wp), DIMENSION(jpij) :: zdivu_adv ! divu as implied by transport scheme (1/s) 137 135 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 138 136 ! 139 INTEGER, PARAMETER :: nitermax = 20137 INTEGER, PARAMETER :: jp_itermax = 20 140 138 !!------------------------------------------------------------------- 141 139 ! controls … … 154 152 ! 0) Identify grid cells with ice 155 153 !-------------------------------- 156 n idx = 0 ; idxice(:) = 0154 npti = 0 ; nptidx(:) = 0 157 155 DO jj = 1, jpj 158 156 DO ji = 1, jpi 159 157 IF ( SUM(a_i(ji,jj,:)) > 0._wp ) THEN 160 n idx = nidx+ 1161 idxice( nidx) = (jj - 1) * jpi + ji158 npti = npti + 1 159 nptidx( npti ) = (jj - 1) * jpi + ji 162 160 ENDIF 163 161 END DO 164 162 END DO 165 163 166 IF( n idx> 0 ) THEN164 IF( npti > 0 ) THEN 167 165 168 166 ! just needed here 169 CALL tab_2d_1d( n idx, idxice(1:nidx), zdivu(1:nidx), divu_i(:,:) )170 CALL tab_2d_1d( n idx, idxice(1:nidx), zdelt(1:nidx), delta_i(:,:) )167 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu(1:npti), divu_i(:,:) ) 168 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt(1:npti), delta_i(:,:) ) 171 169 ! needed here and in the iteration loop 172 CALL tab_3d_2d( n idx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i(:,:,:) )173 CALL tab_3d_2d( n idx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i(:,:,:) )174 CALL tab_2d_1d( n idx, idxice(1:nidx), ato_i_1d(1:nidx) , ato_i(:,:) )175 176 DO ji = 1, n idx170 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i(:,:,:) ) 171 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i(:,:,:) ) 172 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i(:,:) ) 173 174 DO ji = 1, npti 177 175 !-----------------------------------------------------------------------------! 178 176 ! 2) Dynamical inputs (closing rate, divu_adv, opning) … … 219 217 !------------------------------------------------ 220 218 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) 221 nidx2 = 0 ; idxice2(:) = 0222 DO ji = 1, n idx219 ipti = 0 ; iptidx(:) = 0 220 DO ji = 1, npti 223 221 IF( SUM( ABS(apartf(ji,1:jpl)) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 224 nidx2 = nidx2+ 1225 i dxice2 (nidx2) = idxice(ji)226 a_i_2d ( nidx2,:) = a_i_2d (ji,:) ! adjust to new indices227 v_i_2d ( nidx2,:) = v_i_2d (ji,:) ! adjust to new indices228 ato_i_1d ( nidx2) = ato_i_1d (ji) ! adjust to new indices229 closing_net( nidx2) = closing_net(ji) ! adjust to new indices230 zdivu_adv ( nidx2) = zdivu_adv (ji) ! adjust to new indices231 opning ( nidx2) = opning (ji) ! adjust to new indices222 ipti = ipti + 1 223 iptidx (ipti) = nptidx (ji) 224 a_i_2d (ipti,:) = a_i_2d (ji,:) ! adjust to new indices 225 v_i_2d (ipti,:) = v_i_2d (ji,:) ! adjust to new indices 226 ato_i_1d (ipti) = ato_i_1d (ji) ! adjust to new indices 227 closing_net(ipti) = closing_net(ji) ! adjust to new indices 228 zdivu_adv (ipti) = zdivu_adv (ji) ! adjust to new indices 229 opning (ipti) = opning (ji) ! adjust to new indices 232 230 ENDIF 233 231 END DO 234 idxice(:) = idxice2(:)235 n idx = nidx2232 nptidx(:) = iptidx(:) 233 npti = ipti 236 234 237 235 ENDIF … … 240 238 ! Start ridging 241 239 !-------------- 242 IF( n idx> 0 ) THEN240 IF( npti > 0 ) THEN 243 241 244 242 !----------- … … 246 244 !----------- 247 245 ! fields used but not modified 248 CALL tab_2d_1d( n idx, idxice(1:nidx), sss_1d(1:nidx), sss_m(:,:) )249 CALL tab_2d_1d( n idx, idxice(1:nidx), sst_1d(1:nidx), sst_m(:,:) )246 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) ) 247 CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) ) 250 248 ! the following fields are modified in this routine 251 !!CALL tab_2d_1d( n idx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) )252 !!CALL tab_3d_2d( n idx, idxice(1:nidx), a_i_2d(1:nidx,1:jpl), a_i(:,:,:) )253 !!CALL tab_3d_2d( n idx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) )254 CALL tab_3d_2d( n idx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s (:,:,:) )255 CALL tab_3d_2d( n idx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i(:,:,:) )256 CALL tab_3d_2d( n idx, idxice(1:nidx), oa_i_2d(1:nidx,1:jpl), oa_i(:,:,:) )249 !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) ) 250 !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) ) 251 !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 252 CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) ) 253 CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 254 CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 257 255 IF ( nn_pnd_scheme > 0 ) THEN 258 CALL tab_3d_2d( n idx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) )259 CALL tab_3d_2d( n idx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) )256 CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 257 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 260 258 ENDIF 261 259 DO jl = 1, jpl 262 260 DO jk = 1, nlay_s 263 CALL tab_2d_1d( n idx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) )261 CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 264 262 END DO 265 263 DO jk = 1, nlay_i 266 CALL tab_2d_1d( n idx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )267 END DO 268 END DO 269 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_dyn_1d (1:nidx), sfx_dyn (:,:) )270 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri (:,:) )271 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_dyn_1d (1:nidx), wfx_dyn (:,:) )272 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_dyn_1d (1:nidx), hfx_dyn (:,:) )273 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) )264 CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 265 END DO 266 END DO 267 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) ) 268 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) ) 269 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) ) 270 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) ) 271 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 274 272 IF ( nn_pnd_scheme > 0 ) THEN 275 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) )273 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 276 274 ENDIF 277 275 … … 279 277 ! 3) Ridging iteration 280 278 !-----------------------------------------------------------------------------! 281 niter= 1279 iter = 1 282 280 iterate_ridging = 1 283 DO WHILE( iterate_ridging > 0 .AND. niter < nitermax )281 DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax ) 284 282 285 283 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) … … 294 292 ! rates were reduced above), ridge again with new rates. 295 293 iterate_ridging = 0 296 DO ji = 1, n idx294 DO ji = 1, npti 297 295 zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) ) 298 296 IF( ABS( zfac ) < epsi10 ) THEN … … 308 306 END DO 309 307 ! 310 niter = niter + 1311 IF( niter > nitermax ) CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' )308 iter = iter + 1 309 IF( iter > jp_itermax ) CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' ) 312 310 ! 313 311 END DO … … 316 314 ! 1D <==> 2D 317 315 !----------- 318 CALL tab_1d_2d( n idx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) )319 CALL tab_2d_3d( n idx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i (:,:,:) )320 CALL tab_2d_3d( n idx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) )321 CALL tab_2d_3d( n idx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s (:,:,:) )322 CALL tab_2d_3d( n idx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i(:,:,:) )323 CALL tab_2d_3d( n idx, idxice(1:nidx), oa_i_2d(1:nidx,1:jpl), oa_i(:,:,:) )316 CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) ) 317 CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) ) 318 CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 319 CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) ) 320 CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 321 CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 324 322 IF ( nn_pnd_scheme > 0 ) THEN 325 CALL tab_2d_3d( n idx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) )326 CALL tab_2d_3d( n idx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) )323 CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 324 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 327 325 ENDIF 328 326 DO jl = 1, jpl 329 327 DO jk = 1, nlay_s 330 CALL tab_1d_2d( n idx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) )328 CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 331 329 END DO 332 330 DO jk = 1, nlay_i 333 CALL tab_1d_2d( n idx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )334 END DO 335 END DO 336 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_dyn_1d (1:nidx), sfx_dyn (:,:) )337 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri (:,:) )338 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_dyn_1d (1:nidx), wfx_dyn (:,:) )339 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_dyn_1d (1:nidx), hfx_dyn (:,:) )340 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) )331 CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 332 END DO 333 END DO 334 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) ) 335 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) ) 336 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) ) 337 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) ) 338 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 341 339 IF ( nn_pnd_scheme > 0 ) THEN 342 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) )340 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 343 341 ENDIF 344 342 345 ENDIF ! n idx> 0343 ENDIF ! npti > 0 346 344 347 345 CALL ice_var_agg( 1 ) … … 378 376 379 377 ! ! Ice thickness needed for rafting 380 WHERE( pa_i(1:n idx,:) > 0._wp ) ; zhi(1:nidx,:) = pv_i(1:nidx,:) / pa_i(1:nidx,:)381 ELSEWHERE ; zhi(1:n idx,:) = 0._wp378 WHERE( pa_i(1:npti,:) > 0._wp ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 379 ELSEWHERE ; zhi(1:npti,:) = 0._wp 382 380 END WHERE 383 381 … … 398 396 ! Compute total area of ice plus open water. 399 397 ! This is in general not equal to one because of divergence during transport 400 zasum(1:n idx) = pato_i(1:nidx) + SUM( pa_i(1:nidx,:), dim=2 )401 ! 402 WHERE( zasum(1:n idx) > 0._wp ) ; z1_asum(1:nidx) = 1._wp / zasum(1:nidx)403 ELSEWHERE ; z1_asum(1:n idx) = 0._wp398 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 399 ! 400 WHERE( zasum(1:npti) > 0._wp ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 401 ELSEWHERE ; z1_asum(1:npti) = 0._wp 404 402 END WHERE 405 403 ! Compute cumulative thickness distribution function … … 407 405 ! where zGsum(n) is the fractional area in categories 0 to n. 408 406 ! initial value (in h = 0) equals open water area 409 zGsum(1:n idx,-1) = 0._wp410 zGsum(1:n idx,0 ) = pato_i(1:nidx) * z1_asum(1:nidx)407 zGsum(1:npti,-1) = 0._wp 408 zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti) 411 409 DO jl = 1, jpl 412 zGsum(1:n idx,jl) = ( pato_i(1:nidx) + SUM( pa_i(1:nidx,1:jl), dim=2 ) ) * z1_asum(1:nidx) ! sum(1:jl) is ok (and not jpl)410 zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti) ! sum(1:jl) is ok (and not jpl) 413 411 END DO 414 412 ! 415 413 IF( ln_partf_lin ) THEN !--- Linear formulation (Thorndike et al., 1975) 416 414 DO jl = 0, jpl 417 DO ji = 1, n idx415 DO ji = 1, npti 418 416 IF ( zGsum(ji,jl) < rn_gstar ) THEN 419 417 apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * & … … 432 430 zfac = 1._wp / ( 1._wp - EXP(-z1_astar) ) 433 431 DO jl = -1, jpl 434 DO ji = 1, n idx432 DO ji = 1, npti 435 433 zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac 436 434 END DO 437 435 END DO 438 436 DO jl = 0, jpl 439 DO ji = 1, n idx437 DO ji = 1, npti 440 438 apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl) 441 439 END DO … … 447 445 IF( ln_rafting .AND. ln_ridging ) THEN !- ridging & rafting 448 446 DO jl = 1, jpl 449 DO ji = 1, n idx447 DO ji = 1, npti 450 448 aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl) 451 449 araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl) … … 454 452 ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN !- ridging alone 455 453 DO jl = 1, jpl 456 DO ji = 1, n idx454 DO ji = 1, npti 457 455 aridge(ji,jl) = apartf(ji,jl) 458 456 araft (ji,jl) = 0._wp … … 461 459 ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN !- rafting alone 462 460 DO jl = 1, jpl 463 DO ji = 1, n idx461 DO ji = 1, npti 464 462 aridge(ji,jl) = 0._wp 465 463 araft (ji,jl) = apartf(ji,jl) … … 468 466 ELSE !- no ridging & no rafting 469 467 DO jl = 1, jpl 470 DO ji = 1, n idx468 DO ji = 1, npti 471 469 aridge(ji,jl) = 0._wp 472 470 araft (ji,jl) = 0._wp … … 502 500 503 501 zfac = 1._wp / hi_hrft 504 zaksum(1:n idx) = apartf(1:nidx,0)502 zaksum(1:npti) = apartf(1:npti,0) 505 503 ! Transfer function 506 504 DO jl = 1, jpl !all categories have a specific transfer function 507 DO ji = 1, n idx505 DO ji = 1, npti 508 506 IF ( apartf(ji,jl) > 0._wp ) THEN 509 507 zhmean = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min ) … … 530 528 ! closing rate to a gross closing rate. 531 529 ! NOTE: 0 < aksum <= 1 532 WHERE( zaksum(1:n idx) > 0._wp ) ; closing_gross(1:nidx) = closing_net(1:nidx) / zaksum(1:nidx)533 ELSEWHERE ; closing_gross(1:n idx) = 0._wp530 WHERE( zaksum(1:npti) > 0._wp ) ; closing_gross(1:npti) = closing_net(1:npti) / zaksum(1:npti) 531 ELSEWHERE ; closing_gross(1:npti) = 0._wp 534 532 END WHERE 535 533 536 DO ji = 1, n idx534 DO ji = 1, npti 537 535 ! correction to closing rate and opening if closing rate is excessive 538 536 !--------------------------------------------------------------------- … … 552 550 ! would be removed. Reduce the opening rate proportionately. 553 551 DO jl = 1, jpl 554 DO ji = 1, n idx552 DO ji = 1, npti 555 553 zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 556 554 !! IF( zfac > pa_i(ji,jl) .AND. zfac > 0._wp ) THEN … … 594 592 ! 1) Compute change in open water area due to closing and opening. 595 593 !------------------------------------------------------------------------------- 596 DO ji = 1, n idx594 DO ji = 1, npti 597 595 ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice ) 598 596 END DO … … 603 601 DO jl1 = 1, jpl 604 602 605 CALL tab_2d_1d( n idx, idxice(1:nidx), s_i_1d(1:nidx), s_i(:,:,jl1) )606 607 DO ji = 1, n idx603 CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) ) 604 605 DO ji = 1, npti 608 606 609 607 !------------------------------------------------ … … 716 714 ! special loop for e_i because of layers jk 717 715 DO jk = 1, nlay_i 718 DO ji = 1, n idx716 DO ji = 1, npti 719 717 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 720 718 ! Compute ridging /rafting fractions … … 736 734 DO jl2 = 1, jpl 737 735 ! 738 DO ji = 1, n idx736 DO ji = 1, npti 739 737 740 738 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN … … 783 781 784 782 DO jk = 1, nlay_i 785 DO ji = 1, n idx783 DO ji = 1, npti 786 784 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) & 787 785 & ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji) … … 794 792 ! 795 793 ! In case ridging/rafting lead to very small negative values (sometimes it happens) 796 WHERE( a_i_2d(1:n idx,:) < 0._wp ) a_i_2d(1:nidx,:) = 0._wp797 WHERE( v_i_2d(1:n idx,:) < 0._wp ) v_i_2d(1:nidx,:) = 0._wp794 WHERE( a_i_2d(1:npti,:) < 0._wp ) a_i_2d(1:npti,:) = 0._wp 795 WHERE( v_i_2d(1:npti,:) < 0._wp ) v_i_2d(1:npti,:) = 0._wp 798 796 ! 799 797 END SUBROUTINE rdgrft_shift … … 820 818 !!---------------------------------------------------------------------- 821 819 ! !--------------------------------------------------! 822 IF( ln_str_R75 ) THEN ! Ice strength => Rothrock (1975) method ! 823 ! !--------------------------------------------------! 824 ! !--------------------------------------------------! 825 ELSEIF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method ! 820 IF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method ! 826 821 ! !--------------------------------------------------! 827 822 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 828 !829 823 ismooth = 1 830 ! 824 ! !--------------------------------------------------! 825 ELSE ! Zero strength ! 826 ! !--------------------------------------------------! 827 strength(:,:) = 0._wp 828 ismooth = 0 831 829 ENDIF 832 830 ! !--------------------------------------------------! … … 896 894 !! 897 895 NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, & 898 & ln_str_R75, rn_perdg, &899 896 & rn_csrdg , & 900 897 & ln_partf_lin, rn_gstar, & … … 921 918 WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar 922 919 WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg 923 WRITE(numout,*) ' ice strength parameterization Rothrock (1975) ln_str_R75 = ', ln_str_R75924 WRITE(numout,*) ' Ratio of ridging work to PotEner change in ridging rn_perdg = ', rn_perdg925 920 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg 926 921 WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin … … 940 935 ENDIF 941 936 ! 942 IF ( ( ln_str_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN943 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one formulation for ice strength (ln_str_H79 or ln_str_R75)' )944 ENDIF945 !946 937 IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN 947 938 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
r8564 r8565 60 60 ! 61 61 INTEGER :: ji, jj, jl, jcat ! dummy loop index 62 INTEGER :: nidx2! local integer62 INTEGER :: ipti ! local integer 63 63 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 64 64 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - … … 66 66 REAL(wp) :: zslope ! used to compute local thermodynamic "speeds" 67 67 68 INTEGER , DIMENSION(jpij) :: i dxice2! compute remapping or not68 INTEGER , DIMENSION(jpij) :: iptidx ! compute remapping or not 69 69 INTEGER , DIMENSION(jpij,jpl-1) :: jdonor ! donor category index 70 70 REAL(wp), DIMENSION(jpij,jpl) :: zdhice ! ice thickness increment … … 83 83 ! 1) Identify grid cells with ice 84 84 !----------------------------------------------------------------------------------------------- 85 n idx = 0 ; idxice(:) = 085 npti = 0 ; nptidx(:) = 0 86 86 DO jj = 1, jpj 87 87 DO ji = 1, jpi 88 88 IF ( at_i(ji,jj) > epsi10 ) THEN 89 n idx = nidx+ 190 idxice( nidx) = (jj - 1) * jpi + ji89 npti = npti + 1 90 nptidx( npti ) = (jj - 1) * jpi + ji 91 91 ENDIF 92 92 END DO … … 96 96 ! 2) Compute new category boundaries 97 97 !----------------------------------------------------------------------------------------------- 98 IF( n idx> 0 ) THEN98 IF( npti > 0 ) THEN 99 99 100 100 zdhice(:,:) = 0._wp 101 101 zhbnew(:,:) = 0._wp 102 102 103 CALL tab_3d_2d( n idx, idxice(1:nidx), h_i_2d (1:nidx,1:jpl), h_i )104 CALL tab_3d_2d( n idx, idxice(1:nidx), h_ib_2d(1:nidx,1:jpl), h_i_b )105 CALL tab_3d_2d( n idx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i )106 CALL tab_3d_2d( n idx, idxice(1:nidx), a_ib_2d (1:nidx,1:jpl), a_i_b )103 CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i ) 104 CALL tab_3d_2d( npti, nptidx(1:npti), h_ib_2d(1:npti,1:jpl), h_i_b ) 105 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 106 CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d (1:npti,1:jpl), a_i_b ) 107 107 108 108 DO jl = 1, jpl 109 109 ! Compute thickness change in each ice category 110 DO ji = 1, n idx110 DO ji = 1, npti 111 111 zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl) 112 112 END DO … … 116 116 DO jl = 1, jpl - 1 117 117 ! 118 DO ji = 1, n idx118 DO ji = 1, npti 119 119 ! 120 120 ! --- New boundary: Hn* = Hn + Fn*dt --- ! … … 136 136 ! Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 137 137 ! in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 138 IF( a_i_2d(ji,jl ) > epsi10 .AND. h_i_2d(ji,jl ) > ( zhbnew(ji,jl) - epsi10 ) ) idxice(ji) = 0139 IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) ) idxice(ji) = 0138 IF( a_i_2d(ji,jl ) > epsi10 .AND. h_i_2d(ji,jl ) > ( zhbnew(ji,jl) - epsi10 ) ) nptidx(ji) = 0 139 IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) ) nptidx(ji) = 0 140 140 141 141 ! 2) Hn-1 < Hn* < Hn+1 142 IF( zhbnew(ji,jl) < hi_max(jl-1) ) idxice(ji) = 0143 IF( zhbnew(ji,jl) > hi_max(jl+1) ) idxice(ji) = 0142 IF( zhbnew(ji,jl) < hi_max(jl-1) ) nptidx(ji) = 0 143 IF( zhbnew(ji,jl) > hi_max(jl+1) ) nptidx(ji) = 0 144 144 145 145 END DO … … 147 147 ! 148 148 ! --- New boundaries for category jpl --- ! 149 DO ji = 1, n idx149 DO ji = 1, npti 150 150 IF( a_i_2d(ji,jpl) > epsi10 ) THEN 151 151 zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * h_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) ) … … 158 158 ! h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 159 159 ! in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 160 IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) ) idxice(ji) = 0161 IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) ) idxice(ji) = 0160 IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) ) nptidx(ji) = 0 161 IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) ) nptidx(ji) = 0 162 162 END DO 163 163 ! … … 165 165 ! 3) Identify cells where remapping 166 166 !----------------------------------------------------------------------------------------------- 167 nidx2 = 0 ; idxice2(:) = 0168 DO ji = 1, n idx169 IF( idxice(ji) /= 0 ) THEN170 nidx2 = nidx2+ 1171 i dxice2(nidx2) = idxice(ji)172 zhbnew( nidx2,:) = zhbnew(ji,:) ! adjust zhbnew to new indices167 ipti = 0 ; iptidx(:) = 0 168 DO ji = 1, npti 169 IF( nptidx(ji) /= 0 ) THEN 170 ipti = ipti + 1 171 iptidx(ipti) = nptidx(ji) 172 zhbnew(ipti,:) = zhbnew(ji,:) ! adjust zhbnew to new indices 173 173 ENDIF 174 174 END DO 175 idxice(:) = idxice2(:)176 n idx = nidx2175 nptidx(:) = iptidx(:) 176 npti = ipti 177 177 ! 178 178 ENDIF … … 181 181 ! 4) Compute g(h) 182 182 !----------------------------------------------------------------------------------------------- 183 IF( n idx> 0 ) THEN183 IF( npti > 0 ) THEN 184 184 ! 185 185 zhb0(:) = hi_max(0) ; zhb1(:) = hi_max(1) … … 189 189 DO jl = 1, jpl 190 190 ! 191 CALL tab_2d_1d( n idx, idxice(1:nidx), h_ib_1d(1:nidx), h_i_b(:,:,jl) )192 CALL tab_2d_1d( n idx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,jl) )193 CALL tab_2d_1d( n idx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) )194 CALL tab_2d_1d( n idx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) )191 CALL tab_2d_1d( npti, nptidx(1:npti), h_ib_1d(1:npti), h_i_b(:,:,jl) ) 192 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl) ) 193 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i(:,:,jl) ) 194 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i(:,:,jl) ) 195 195 ! 196 196 IF( jl == 1 ) THEN 197 197 ! 198 198 ! --- g(h) for category 1 --- ! 199 CALL ice_itd_glinear( zhb0(1:n idx) , zhb1(1:nidx) , h_ib_1d(1:nidx) , a_i_1d(1:nidx) , & ! in200 & g0 (1:n idx,1), g1 (1:nidx,1), hL (1:nidx,1), hR (1:nidx,1) ) ! out199 CALL ice_itd_glinear( zhb0(1:npti) , zhb1(1:npti) , h_ib_1d(1:npti) , a_i_1d(1:npti) , & ! in 200 & g0 (1:npti,1), g1 (1:npti,1), hL (1:npti,1), hR (1:npti,1) ) ! out 201 201 ! 202 202 ! Area lost due to melting of thin ice 203 DO ji = 1, n idx203 DO ji = 1, npti 204 204 ! 205 205 IF( a_i_1d(ji) > epsi10 ) THEN … … 233 233 END DO 234 234 ! 235 CALL tab_1d_2d( n idx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,jl) )236 CALL tab_1d_2d( n idx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) )237 CALL tab_1d_2d( n idx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) )235 CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl) ) 236 CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i(:,:,jl) ) 237 CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i(:,:,jl) ) 238 238 ! 239 239 ENDIF ! jl=1 240 240 ! 241 241 ! --- g(h) for each thickness category --- ! 242 CALL ice_itd_glinear( zhbnew(1:n idx,jl-1), zhbnew(1:nidx,jl), h_i_1d(1:nidx) , a_i_1d(1:nidx) , & ! in243 & g0 (1:n idx,jl ), g1 (1:nidx,jl), hL (1:nidx,jl), hR (1:nidx,jl) ) ! out242 CALL ice_itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti) , a_i_1d(1:npti) , & ! in 243 & g0 (1:npti,jl ), g1 (1:npti,jl), hL (1:npti,jl), hR (1:npti,jl) ) ! out 244 244 ! 245 245 END DO … … 250 250 DO jl = 1, jpl - 1 251 251 ! 252 DO ji = 1, n idx252 DO ji = 1, npti 253 253 ! 254 254 ! left and right integration limits in eta space … … 281 281 ! 6) Shift ice between categories 282 282 !---------------------------------------------------------------------------------------------- 283 CALL ice_itd_shiftice ( jdonor(1:n idx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )283 CALL ice_itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) 284 284 285 285 !---------------------------------------------------------------------------------------------- 286 286 ! 7) Make sure h_i >= minimum ice thickness hi_min 287 287 !---------------------------------------------------------------------------------------------- 288 CALL tab_2d_1d( n idx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,1) )289 CALL tab_2d_1d( n idx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,1) )290 CALL tab_2d_1d( n idx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1) )288 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,1) ) 289 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i(:,:,1) ) 290 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1) ) 291 291 292 DO ji = 1, n idx292 DO ji = 1, npti 293 293 IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 294 294 a_i_1d (ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin … … 302 302 END DO 303 303 ! 304 CALL tab_1d_2d( n idx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,1) )305 CALL tab_1d_2d( n idx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,1) )306 CALL tab_1d_2d( n idx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1) )304 CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,1) ) 305 CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i(:,:,1) ) 306 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1) ) 307 307 ! 308 308 ENDIF … … 338 338 z2_3 = 2._wp / 3._wp 339 339 ! 340 DO ji = 1, n idx340 DO ji = 1, npti 341 341 ! 342 342 IF( paice(ji) > epsi10 .AND. phice(ji) > 0._wp ) THEN … … 392 392 !!------------------------------------------------------------------ 393 393 394 CALL tab_3d_2d( n idx, idxice(1:nidx), h_i_2d (1:nidx,1:jpl), h_i )395 CALL tab_3d_2d( n idx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i )396 CALL tab_3d_2d( n idx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i )397 CALL tab_3d_2d( n idx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s )398 CALL tab_3d_2d( n idx, idxice(1:nidx), oa_i_2d(1:nidx,1:jpl), oa_i )399 CALL tab_3d_2d( n idx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i )400 CALL tab_3d_2d( n idx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip )401 CALL tab_3d_2d( n idx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip )402 CALL tab_3d_2d( n idx, idxice(1:nidx), t_su_2d(1:nidx,1:jpl), t_su )394 CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i ) 395 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 396 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) 397 CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s ) 398 CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i ) 399 CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i ) 400 CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip ) 401 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 402 CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 403 403 404 404 !---------------------------------------------------------------------------------------------- … … 406 406 !---------------------------------------------------------------------------------------------- 407 407 DO jl = 1, jpl 408 DO ji = 1, n idx408 DO ji = 1, npti 409 409 zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl) 410 410 END DO … … 415 415 !------------------------------------------------------------------------------- 416 416 DO jl = 1, jpl - 1 417 DO ji = 1, n idx417 DO ji = 1, npti 418 418 ! 419 419 jl1 = kdonor(ji,jl) … … 475 475 DO jk = 1, nlay_s !--- Snow heat content 476 476 ! 477 DO ji = 1, n idx478 ii = MOD( idxice(ji) - 1, jpi ) + 1479 ij = ( idxice(ji) - 1 ) / jpi + 1477 DO ji = 1, npti 478 ii = MOD( nptidx(ji) - 1, jpi ) + 1 479 ij = ( nptidx(ji) - 1 ) / jpi + 1 480 480 ! 481 481 jl1 = kdonor(ji,jl) … … 494 494 495 495 DO jk = 1, nlay_i !--- Ice heat content 496 DO ji = 1, n idx497 ii = MOD( idxice(ji) - 1, jpi ) + 1498 ij = ( idxice(ji) - 1 ) / jpi + 1496 DO ji = 1, npti 497 ii = MOD( nptidx(ji) - 1, jpi ) + 1 498 ij = ( nptidx(ji) - 1 ) / jpi + 1 499 499 ! 500 500 jl1 = kdonor(ji,jl) … … 517 517 ! 3) Update ice thickness and temperature 518 518 !------------------------------------------------------------------------------- 519 WHERE( a_i_2d(1:n idx,:) >= epsi20 )520 h_i_2d(1:n idx,:) = v_i_2d(1:nidx,:) / a_i_2d(1:nidx,:)521 t_su_2d(1:n idx,:) = zaTsfn(1:nidx,:) / a_i_2d(1:nidx,:)519 WHERE( a_i_2d(1:npti,:) >= epsi20 ) 520 h_i_2d(1:npti,:) = v_i_2d(1:npti,:) / a_i_2d(1:npti,:) 521 t_su_2d(1:npti,:) = zaTsfn(1:npti,:) / a_i_2d(1:npti,:) 522 522 ELSEWHERE 523 h_i_2d(1:n idx,:) = 0._wp524 t_su_2d(1:n idx,:) = rt0523 h_i_2d(1:npti,:) = 0._wp 524 t_su_2d(1:npti,:) = rt0 525 525 END WHERE 526 526 ! 527 CALL tab_2d_3d( n idx, idxice(1:nidx), h_i_2d (1:nidx,1:jpl), h_i )528 CALL tab_2d_3d( n idx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i )529 CALL tab_2d_3d( n idx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i )530 CALL tab_2d_3d( n idx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s )531 CALL tab_2d_3d( n idx, idxice(1:nidx), oa_i_2d(1:nidx,1:jpl), oa_i )532 CALL tab_2d_3d( n idx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i )533 CALL tab_2d_3d( n idx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip )534 CALL tab_2d_3d( n idx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip )535 CALL tab_2d_3d( n idx, idxice(1:nidx), t_su_2d(1:nidx,1:jpl), t_su )527 CALL tab_2d_3d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i ) 528 CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 529 CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) 530 CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s ) 531 CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i ) 532 CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i ) 533 CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip ) 534 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 535 CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 536 536 ! 537 537 END SUBROUTINE ice_itd_shiftice … … 564 564 DO jl = 1, jpl-1 ! identify thicknesses that are too big 565 565 ! !--------------------------------------- 566 n idx = 0 ; idxice(:) = 0566 npti = 0 ; nptidx(:) = 0 567 567 DO jj = 1, jpj 568 568 DO ji = 1, jpi 569 569 IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN 570 n idx = nidx+ 1571 idxice( nidx) = (jj - 1) * jpi + ji570 npti = npti + 1 571 nptidx( npti ) = (jj - 1) * jpi + ji 572 572 ENDIF 573 573 END DO 574 574 END DO 575 575 ! 576 !!clem CALL tab_2d_1d( n idx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,jl) )577 CALL tab_2d_1d( n idx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) )578 CALL tab_2d_1d( n idx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) )579 ! 580 DO ji = 1, n idx576 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl) ) 577 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i(:,:,jl) ) 578 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i(:,:,jl) ) 579 ! 580 DO ji = 1, npti 581 581 jdonor(ji,jl) = jl 582 582 ! how much of a_i you send in cat sup is somewhat arbitrary … … 592 592 END DO 593 593 ! 594 IF( n idx> 0 ) THEN595 CALL ice_itd_shiftice( jdonor(1:n idx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) ) ! Shift jl=>jl+1594 IF( npti > 0 ) THEN 595 CALL ice_itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl=>jl+1 596 596 ! Reset shift parameters 597 jdonor(1:n idx,jl) = 0598 zdaice(1:n idx,jl) = 0._wp599 zdvice(1:n idx,jl) = 0._wp597 jdonor(1:npti,jl) = 0 598 zdaice(1:npti,jl) = 0._wp 599 zdvice(1:npti,jl) = 0._wp 600 600 ENDIF 601 601 ! … … 605 605 DO jl = jpl-1, 1, -1 ! Identify thicknesses that are too small 606 606 ! !----------------------------------------- 607 n idx = 0 ; idxice(:) = 0607 npti = 0 ; nptidx(:) = 0 608 608 DO jj = 1, jpj 609 609 DO ji = 1, jpi 610 610 IF( a_i(ji,jj,jl+1) > epsi10 .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN 611 n idx = nidx+ 1612 idxice( nidx) = (jj - 1) * jpi + ji611 npti = npti + 1 612 nptidx( npti ) = (jj - 1) * jpi + ji 613 613 ENDIF 614 614 END DO 615 615 END DO 616 616 ! 617 CALL tab_2d_1d( n idx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl+1) ) ! jl+1 is ok618 CALL tab_2d_1d( n idx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl+1) ) ! jl+1 is ok619 DO ji = 1, n idx617 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 618 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 619 DO ji = 1, npti 620 620 jdonor(ji,jl) = jl + 1 621 621 zdaice(ji,jl) = a_i_1d(ji) … … 623 623 END DO 624 624 ! 625 IF( n idx> 0 ) THEN626 CALL ice_itd_shiftice( jdonor(1:n idx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) ) ! Shift jl+1=>jl625 IF( npti > 0 ) THEN 626 CALL ice_itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl+1=>jl 627 627 ! Reset shift parameters 628 jdonor(1:n idx,jl) = 0629 zdaice(1:n idx,jl) = 0._wp630 zdvice(1:n idx,jl) = 0._wp628 jdonor(1:npti,jl) = 0 629 zdaice(1:npti,jl) = 0._wp 630 zdvice(1:npti,jl) = 0._wp 631 631 ENDIF 632 632 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
r8564 r8565 211 211 212 212 ! select ice covered grid points 213 n idx = 0 ; idxice(:) = 0213 npti = 0 ; nptidx(:) = 0 214 214 DO jj = 1, jpj 215 215 DO ji = 1, jpi 216 216 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 217 n idx = nidx+ 1218 idxice(nidx) = (jj - 1) * jpi + ji217 npti = npti + 1 218 nptidx(npti) = (jj - 1) * jpi + ji 219 219 ENDIF 220 220 END DO 221 221 END DO 222 222 223 IF( lk_mpp ) CALL mpp_ini_ice( n idx, numout )224 225 IF( n idx> 0 ) THEN ! If there is no ice, do nothing.223 IF( lk_mpp ) CALL mpp_ini_ice( npti , numout ) 224 225 IF( npti > 0 ) THEN ! If there is no ice, do nothing. 226 226 ! 227 227 CALL ice_thd_1d2d( jl, 1 ) ! --- Move to 1D arrays --- ! 228 228 ! ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- ! 229 229 ! 230 s_i_new (1:n idx) = 0._wp ; dh_s_tot (1:nidx) = 0._wp ! --- some init --- ! (important to have them here)231 dh_i_surf (1:n idx) = 0._wp ; dh_i_bott(1:nidx) = 0._wp232 dh_snowice(1:n idx) = 0._wp ; dh_i_sub (1:nidx) = 0._wp230 s_i_new (1:npti) = 0._wp ; dh_s_tot (1:npti) = 0._wp ! --- some init --- ! (important to have them here) 231 dh_i_surf (1:npti) = 0._wp ; dh_i_bott(1:npti) = 0._wp 232 dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp 233 233 ! 234 234 IF( ln_icedH ) THEN ! --- growing/melting --- ! 235 235 CALL ice_thd_zdf ! Ice/Snow Temperature profile 236 236 CALL ice_thd_dh ! Ice/Snow thickness 237 CALL ice_thd_ent( e_i_1d(1:n idx,:) ) ! Ice enthalpy remapping237 CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping 238 238 ENDIF 239 239 ! … … 294 294 ! Recover ice temperature 295 295 DO jk = 1, nlay_i 296 DO ji = 1, n idx296 DO ji = 1, npti 297 297 ztmelts = -tmut * sz_i_1d(ji,jk) 298 298 ! Conversion q(S,T) -> T (second order equation) … … 323 323 !!----------------------------------------------------------------------- 324 324 ! 325 DO ji = 1, n idx325 DO ji = 1, npti 326 326 zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) + dh_i_sub(ji) ) 327 327 IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp ) THEN … … 360 360 CASE( 1 ) !== from 2D to 1D ==! 361 361 ! !---------------------! 362 CALL tab_2d_1d( n idx, idxice(1:nidx), at_i_1d(1:nidx), at_i )363 CALL tab_2d_1d( n idx, idxice(1:nidx), a_i_1d (1:nidx), a_i (:,:,kl) )364 CALL tab_2d_1d( n idx, idxice(1:nidx), h_i_1d(1:nidx), h_i(:,:,kl) )365 CALL tab_2d_1d( n idx, idxice(1:nidx), h_s_1d(1:nidx), h_s(:,:,kl) )366 CALL tab_2d_1d( n idx, idxice(1:nidx), t_su_1d(1:nidx), t_su(:,:,kl) )367 CALL tab_2d_1d( n idx, idxice(1:nidx), s_i_1d(1:nidx), s_i(:,:,kl) )362 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i ) 363 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl) ) 364 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl) ) 365 CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl) ) 366 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl) ) 367 CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl) ) 368 368 DO jk = 1, nlay_s 369 CALL tab_2d_1d( n idx, idxice(1:nidx), t_s_1d(1:nidx,jk), t_s(:,:,jk,kl) )370 CALL tab_2d_1d( n idx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,kl) )369 CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) ) 370 CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) ) 371 371 END DO 372 372 DO jk = 1, nlay_i 373 CALL tab_2d_1d( n idx, idxice(1:nidx), t_i_1d(1:nidx,jk), t_i(:,:,jk,kl) )374 CALL tab_2d_1d( n idx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,kl) )375 CALL tab_2d_1d( n idx, idxice(1:nidx), sz_i_1d(1:nidx,jk), sz_i(:,:,jk,kl) )376 END DO 377 ! 378 CALL tab_2d_1d( n idx, idxice(1:nidx), qprec_ice_1d(1:nidx), qprec_ice )379 CALL tab_2d_1d( n idx, idxice(1:nidx), qsr_ice_1d (1:nidx), qsr_ice (:,:,kl) )380 CALL tab_2d_1d( n idx, idxice(1:nidx), fr1_i0_1d (1:nidx), fr1_i0 )381 CALL tab_2d_1d( n idx, idxice(1:nidx), fr2_i0_1d (1:nidx), fr2_i0 )382 CALL tab_2d_1d( n idx, idxice(1:nidx), qns_ice_1d (1:nidx), qns_ice (:,:,kl) )383 CALL tab_2d_1d( n idx, idxice(1:nidx), ftr_ice_1d (1:nidx), ftr_ice (:,:,kl) )384 CALL tab_2d_1d( n idx, idxice(1:nidx), evap_ice_1d (1:nidx), evap_ice(:,:,kl) )385 CALL tab_2d_1d( n idx, idxice(1:nidx), dqns_ice_1d (1:nidx), dqns_ice(:,:,kl) )386 CALL tab_2d_1d( n idx, idxice(1:nidx), t_bo_1d (1:nidx), t_bo )387 CALL tab_2d_1d( n idx, idxice(1:nidx), sprecip_1d (1:nidx), sprecip )388 CALL tab_2d_1d( n idx, idxice(1:nidx), fhtur_1d (1:nidx), fhtur )389 CALL tab_2d_1d( n idx, idxice(1:nidx), fhld_1d (1:nidx), fhld )390 ! 391 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_snw_sni_1d(1:nidx), wfx_snw_sni )392 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_snw_sum_1d(1:nidx), wfx_snw_sum )393 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_sub_1d (1:nidx), wfx_sub )394 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_snw_sub_1d(1:nidx), wfx_snw_sub )395 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_ice_sub_1d(1:nidx), wfx_ice_sub )396 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_err_sub_1d(1:nidx), wfx_err_sub )397 ! 398 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_bog_1d (1:nidx), wfx_bog )399 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_bom_1d (1:nidx), wfx_bom )400 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_sum_1d (1:nidx), wfx_sum )401 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_sni_1d (1:nidx), wfx_sni )402 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_res_1d (1:nidx), wfx_res )403 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_spr_1d (1:nidx), wfx_spr )404 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_lam_1d (1:nidx), wfx_lam )405 ! 406 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_bog_1d (1:nidx), sfx_bog )407 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_bom_1d (1:nidx), sfx_bom )408 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_sum_1d (1:nidx), sfx_sum )409 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_sni_1d (1:nidx), sfx_sni )410 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri )411 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_res_1d (1:nidx), sfx_res )412 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_sub_1d (1:nidx), sfx_sub )413 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_lam_1d (1:nidx), sfx_lam )414 ! 415 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_thd_1d (1:nidx), hfx_thd )416 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_spr_1d (1:nidx), hfx_spr )417 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_sum_1d (1:nidx), hfx_sum )418 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_bom_1d (1:nidx), hfx_bom )419 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_bog_1d (1:nidx), hfx_bog )420 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_dif_1d (1:nidx), hfx_dif )421 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_opw_1d (1:nidx), hfx_opw )422 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw )423 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub )424 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res )425 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif )426 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_err_rem_1d(1:nidx), hfx_err_rem )427 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_out_1d (1:nidx), hfx_out )373 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl) ) 374 CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl) ) 375 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 376 END DO 377 ! 378 CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d(1:npti), qprec_ice ) 379 CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_1d (1:npti), qsr_ice (:,:,kl) ) 380 CALL tab_2d_1d( npti, nptidx(1:npti), fr1_i0_1d (1:npti), fr1_i0 ) 381 CALL tab_2d_1d( npti, nptidx(1:npti), fr2_i0_1d (1:npti), fr2_i0 ) 382 CALL tab_2d_1d( npti, nptidx(1:npti), qns_ice_1d (1:npti), qns_ice (:,:,kl) ) 383 CALL tab_2d_1d( npti, nptidx(1:npti), ftr_ice_1d (1:npti), ftr_ice (:,:,kl) ) 384 CALL tab_2d_1d( npti, nptidx(1:npti), evap_ice_1d (1:npti), evap_ice(:,:,kl) ) 385 CALL tab_2d_1d( npti, nptidx(1:npti), dqns_ice_1d (1:npti), dqns_ice(:,:,kl) ) 386 CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d (1:npti), t_bo ) 387 CALL tab_2d_1d( npti, nptidx(1:npti), sprecip_1d (1:npti), sprecip ) 388 CALL tab_2d_1d( npti, nptidx(1:npti), fhtur_1d (1:npti), fhtur ) 389 CALL tab_2d_1d( npti, nptidx(1:npti), fhld_1d (1:npti), fhld ) 390 ! 391 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) 392 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 393 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sub_1d (1:npti), wfx_sub ) 394 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub ) 395 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub ) 396 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub ) 397 ! 398 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog ) 399 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom ) 400 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum ) 401 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni ) 402 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res ) 403 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr ) 404 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam ) 405 ! 406 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog ) 407 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom ) 408 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum ) 409 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni ) 410 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri ) 411 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res ) 412 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub ) 413 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam ) 414 ! 415 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d (1:npti), hfx_thd ) 416 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_spr_1d (1:npti), hfx_spr ) 417 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sum_1d (1:npti), hfx_sum ) 418 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bom_1d (1:npti), hfx_bom ) 419 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bog_1d (1:npti), hfx_bog ) 420 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dif_1d (1:npti), hfx_dif ) 421 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d (1:npti), hfx_opw ) 422 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_snw_1d (1:npti), hfx_snw ) 423 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sub_1d (1:npti), hfx_sub ) 424 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 425 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 426 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem ) 427 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_out_1d (1:npti), hfx_out ) 428 428 ! 429 429 ! SIMIP diagnostics 430 CALL tab_2d_1d( n idx, idxice(1:nidx), diag_fc_bo_1d(1:nidx), diag_fc_bo )431 CALL tab_2d_1d( n idx, idxice(1:nidx), diag_fc_su_1d(1:nidx), diag_fc_su )430 CALL tab_2d_1d( npti, nptidx(1:npti), diag_fc_bo_1d(1:npti), diag_fc_bo ) 431 CALL tab_2d_1d( npti, nptidx(1:npti), diag_fc_su_1d(1:npti), diag_fc_su ) 432 432 ! ocean surface fields 433 CALL tab_2d_1d( n idx, idxice(1:nidx), sst_1d(1:nidx), sst_m )434 CALL tab_2d_1d( n idx, idxice(1:nidx), sss_1d(1:nidx), sss_m )433 CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m ) 434 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m ) 435 435 436 436 ! --- Change units of e_i, e_s from J/m2 to J/m3 --- ! 437 437 DO jk = 1, nlay_i 438 WHERE( h_i_1d(1:n idx)>0._wp ) e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) / (h_i_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_i438 WHERE( h_i_1d(1:npti)>0._wp ) e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i 439 439 END DO 440 440 DO jk = 1, nlay_s 441 WHERE( h_s_1d(1:n idx)>0._wp ) e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) / (h_s_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_s441 WHERE( h_s_1d(1:npti)>0._wp ) e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s 442 442 END DO 443 443 ! … … 447 447 ! --- Change units of e_i, e_s from J/m3 to J/m2 --- ! 448 448 DO jk = 1, nlay_i 449 e_i_1d(1:n idx,jk) = e_i_1d(1:nidx,jk) * h_i_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_i449 e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) * h_i_1d(1:npti) * a_i_1d(1:npti) * r1_nlay_i 450 450 END DO 451 451 DO jk = 1, nlay_s 452 e_s_1d(1:n idx,jk) = e_s_1d(1:nidx,jk) * h_s_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_s452 e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) * h_s_1d(1:npti) * a_i_1d(1:npti) * r1_nlay_s 453 453 END DO 454 454 ! 455 455 ! Change thickness to volume (replaces routine ice_var_eqv2glo) 456 v_i_1d(1:n idx) = h_i_1d(1:nidx) * a_i_1d(1:nidx)457 v_s_1d(1:n idx) = h_s_1d(1:nidx) * a_i_1d(1:nidx)458 sv_i_1d(1:n idx) = s_i_1d(1:nidx) * v_i_1d(1:nidx)456 v_i_1d(1:npti) = h_i_1d(1:npti) * a_i_1d(1:npti) 457 v_s_1d(1:npti) = h_s_1d(1:npti) * a_i_1d(1:npti) 458 sv_i_1d(1:npti) = s_i_1d(1:npti) * v_i_1d(1:npti) 459 459 460 CALL tab_1d_2d( n idx, idxice(1:nidx), at_i_1d(1:nidx), at_i )461 CALL tab_1d_2d( n idx, idxice(1:nidx), a_i_1d (1:nidx), a_i (:,:,kl) )462 CALL tab_1d_2d( n idx, idxice(1:nidx), h_i_1d(1:nidx), h_i(:,:,kl) )463 CALL tab_1d_2d( n idx, idxice(1:nidx), h_s_1d(1:nidx), h_s(:,:,kl) )464 CALL tab_1d_2d( n idx, idxice(1:nidx), t_su_1d(1:nidx), t_su(:,:,kl) )465 CALL tab_1d_2d( n idx, idxice(1:nidx), s_i_1d(1:nidx), s_i(:,:,kl) )460 CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i ) 461 CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl) ) 462 CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl) ) 463 CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl) ) 464 CALL tab_1d_2d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl) ) 465 CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl) ) 466 466 DO jk = 1, nlay_s 467 CALL tab_1d_2d( n idx, idxice(1:nidx), t_s_1d(1:nidx,jk), t_s(:,:,jk,kl) )468 CALL tab_1d_2d( n idx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,kl) )467 CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) ) 468 CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) ) 469 469 END DO 470 470 DO jk = 1, nlay_i 471 CALL tab_1d_2d( n idx, idxice(1:nidx), t_i_1d(1:nidx,jk), t_i(:,:,jk,kl) )472 CALL tab_1d_2d( n idx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,kl) )473 CALL tab_1d_2d( n idx, idxice(1:nidx), sz_i_1d(1:nidx,jk), sz_i(:,:,jk,kl) )474 END DO 475 ! 476 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_snw_sni_1d(1:nidx), wfx_snw_sni )477 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_snw_sum_1d(1:nidx), wfx_snw_sum )478 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_sub_1d (1:nidx), wfx_sub )479 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_snw_sub_1d(1:nidx), wfx_snw_sub )480 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_ice_sub_1d(1:nidx), wfx_ice_sub )481 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_err_sub_1d(1:nidx), wfx_err_sub )482 ! 483 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_bog_1d (1:nidx), wfx_bog )484 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_bom_1d (1:nidx), wfx_bom )485 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_sum_1d (1:nidx), wfx_sum )486 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_sni_1d (1:nidx), wfx_sni )487 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_res_1d (1:nidx), wfx_res )488 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_spr_1d (1:nidx), wfx_spr )489 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_lam_1d (1:nidx), wfx_lam )490 ! 491 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_bog_1d (1:nidx), sfx_bog )492 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_bom_1d (1:nidx), sfx_bom )493 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_sum_1d (1:nidx), sfx_sum )494 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_sni_1d (1:nidx), sfx_sni )495 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri )496 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_res_1d (1:nidx), sfx_res )497 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_sub_1d (1:nidx), sfx_sub )498 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_lam_1d (1:nidx), sfx_lam )499 ! 500 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_thd_1d (1:nidx), hfx_thd )501 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_spr_1d (1:nidx), hfx_spr )502 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_sum_1d (1:nidx), hfx_sum )503 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_bom_1d (1:nidx), hfx_bom )504 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_bog_1d (1:nidx), hfx_bog )505 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_dif_1d (1:nidx), hfx_dif )506 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_opw_1d (1:nidx), hfx_opw )507 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw )508 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub )509 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res )510 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif )511 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_err_rem_1d(1:nidx), hfx_err_rem )512 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_out_1d (1:nidx), hfx_out )513 ! 514 CALL tab_1d_2d( n idx, idxice(1:nidx), qns_ice_1d(1:nidx), qns_ice(:,:,kl) )515 CALL tab_1d_2d( n idx, idxice(1:nidx), ftr_ice_1d(1:nidx), ftr_ice(:,:,kl) )471 CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl) ) 472 CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl) ) 473 CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 474 END DO 475 ! 476 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) 477 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 478 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sub_1d (1:npti), wfx_sub ) 479 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub ) 480 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub ) 481 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub ) 482 ! 483 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog ) 484 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom ) 485 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum ) 486 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni ) 487 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res ) 488 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr ) 489 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam ) 490 ! 491 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog ) 492 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom ) 493 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum ) 494 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni ) 495 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri ) 496 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res ) 497 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub ) 498 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam ) 499 ! 500 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d (1:npti), hfx_thd ) 501 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_spr_1d (1:npti), hfx_spr ) 502 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sum_1d (1:npti), hfx_sum ) 503 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bom_1d (1:npti), hfx_bom ) 504 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bog_1d (1:npti), hfx_bog ) 505 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dif_1d (1:npti), hfx_dif ) 506 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d (1:npti), hfx_opw ) 507 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_snw_1d (1:npti), hfx_snw ) 508 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sub_1d (1:npti), hfx_sub ) 509 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 510 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 511 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem ) 512 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_out_1d (1:npti), hfx_out ) 513 ! 514 CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d(1:npti), qns_ice(:,:,kl) ) 515 CALL tab_1d_2d( npti, nptidx(1:npti), ftr_ice_1d(1:npti), ftr_ice(:,:,kl) ) 516 516 ! 517 517 ! SIMIP diagnostics 518 CALL tab_1d_2d( n idx, idxice(1:nidx), t_si_1d (1:nidx), t_si(:,:,kl) )519 CALL tab_1d_2d( n idx, idxice(1:nidx), diag_fc_bo_1d(1:nidx), diag_fc_bo )520 CALL tab_1d_2d( n idx, idxice(1:nidx), diag_fc_su_1d(1:nidx), diag_fc_su )518 CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d (1:npti), t_si(:,:,kl) ) 519 CALL tab_1d_2d( npti, nptidx(1:npti), diag_fc_bo_1d(1:npti), diag_fc_bo ) 520 CALL tab_1d_2d( npti, nptidx(1:npti), diag_fc_su_1d(1:npti), diag_fc_su ) 521 521 ! extensive variables 522 CALL tab_1d_2d( n idx, idxice(1:nidx), v_i_1d (1:nidx), v_i (:,:,kl) )523 CALL tab_1d_2d( n idx, idxice(1:nidx), v_s_1d (1:nidx), v_s (:,:,kl) )524 CALL tab_1d_2d( n idx, idxice(1:nidx), sv_i_1d(1:nidx), sv_i(:,:,kl) )522 CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,kl) ) 523 CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) ) 524 CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) ) 525 525 ! 526 526 END SELECT -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90
r8564 r8565 118 118 ! 119 119 zastar = 1._wp / ( 1._wp - (rn_dmin / zdmax)**(1._wp/rn_beta) ) 120 DO ji = 1, n idx120 DO ji = 1, npti 121 121 ! --- Calculate reduction of total sea ice concentration --- ! 122 122 zdfloe = rn_dmin * ( zastar / ( zastar - at_i_1d(ji) ) )**rn_beta ! Mean floe caliper diameter [m] -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
r8564 r8565 117 117 118 118 ! initialize layer thicknesses and enthalpies 119 h_i_old (1:n idx,0:nlay_i+1) = 0._wp120 eh_i_old(1:n idx,0:nlay_i+1) = 0._wp119 h_i_old (1:npti,0:nlay_i+1) = 0._wp 120 eh_i_old(1:npti,0:nlay_i+1) = 0._wp 121 121 DO jk = 1, nlay_i 122 DO ji = 1, n idx122 DO ji = 1, npti 123 123 h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i 124 124 eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk) … … 130 130 !------------------------------------------------------------------------------! 131 131 ! 132 DO ji = 1, n idx132 DO ji = 1, npti 133 133 zdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 134 134 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) … … 142 142 ! (should not happen but sometimes it does) 143 143 !------------------------------------------------------------------------------! 144 DO ji = 1, n idx144 DO ji = 1, npti 145 145 IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting 146 146 ! Contribution to heat flux to the ocean [W.m-2], < 0 … … 159 159 !------------------------------------------------------------! 160 160 DO jk = 1, nlay_i 161 DO ji = 1, n idx161 DO ji = 1, npti 162 162 zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i 163 163 zeh_i(ji) = zeh_i(ji) + e_i_1d(ji,jk) * zh_i(ji,jk) … … 183 183 ! Martin Vancoppenolle, December 2006 184 184 185 CALL ice_thd_snwblow( 1. - at_i_1d(1:n idx), zsnw(1:nidx) ) ! snow distribution over ice after wind blowing186 187 zdeltah(1:n idx,:) = 0._wp188 DO ji = 1, n idx185 CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing 186 187 zdeltah(1:npti,:) = 0._wp 188 DO ji = 1, npti 189 189 !----------- 190 190 ! Snow fall … … 220 220 221 221 ! If heat still available (zq_su > 0), then melt more snow 222 zdeltah(1:n idx,:) = 0._wp223 zdh_s_mel(1:n idx) = 0._wp222 zdeltah(1:npti,:) = 0._wp 223 zdh_s_mel(1:npti) = 0._wp 224 224 DO jk = 1, nlay_s 225 DO ji = 1, n idx225 DO ji = 1, npti 226 226 ! thickness change 227 227 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) … … 245 245 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 246 246 ! clem comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean) 247 zdeltah(1:n idx,:) = 0._wp248 DO ji = 1, n idx247 zdeltah(1:npti,:) = 0._wp 248 DO ji = 1, npti 249 249 zdh_s_sub(ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 250 250 ! remaining evap in kg.m-2 (used for ice melting later on) … … 265 265 266 266 ! --- Update snow diags --- ! 267 DO ji = 1, n idx267 DO ji = 1, npti 268 268 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 269 269 END DO … … 274 274 ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 275 275 DO jk = 1, nlay_s 276 DO ji = 1,n idx276 DO ji = 1,npti 277 277 rswitch = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) ) 278 278 e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) * & … … 285 285 ! 3.4 Surface ice ablation 286 286 !-------------------------- 287 zdeltah(1:n idx,:) = 0._wp ! important287 zdeltah(1:npti,:) = 0._wp ! important 288 288 DO jk = 1, nlay_i 289 DO ji = 1, n idx289 DO ji = 1, npti 290 290 ztmelts = - tmut * sz_i_1d(ji,jk) ! Melting point of layer k [C] 291 291 … … 375 375 END DO 376 376 ! update ice thickness 377 DO ji = 1, n idx377 DO ji = 1, npti 378 378 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) ) 379 379 END DO 380 380 381 381 ! remaining "potential" evap is sent to ocean 382 DO ji = 1, n idx382 DO ji = 1, npti 383 383 wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice ! <=0 (net evap for the ocean in kg.m-2.s-1) 384 384 END DO … … 407 407 408 408 ! Iterative procedure 409 DO ji = 1, n idx409 DO ji = 1, npti 410 410 IF( zf_tt(ji) < 0._wp ) THEN 411 411 DO iter = 1, num_iter_max … … 478 478 ! 4.2 Basal melt 479 479 !---------------- 480 zdeltah(1:n idx,:) = 0._wp ! important480 zdeltah(1:npti,:) = 0._wp ! important 481 481 DO jk = nlay_i, 1, -1 482 DO ji = 1, n idx482 DO ji = 1, npti 483 483 IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji,jk) ) THEN ! do not calculate where layer has already disappeared by surface melting 484 484 … … 554 554 ! Update temperature, energy 555 555 !------------------------------------------- 556 DO ji = 1, n idx556 DO ji = 1, npti 557 557 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bott(ji) ) 558 558 END DO … … 563 563 ! If heat still available for melting and snow remains, then melt more snow 564 564 !------------------------------------------- 565 zdeltah(1:n idx,:) = 0._wp ! important566 DO ji = 1, n idx565 zdeltah(1:npti,:) = 0._wp ! important 566 DO ji = 1, npti 567 567 zq_rema(ji) = zq_su(ji) + zq_bo(ji) 568 568 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ! =1 if snow … … 592 592 ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 593 593 z1_rho = 1._wp / ( rhosn+rau0-rhoic ) 594 DO ji = 1, n idx594 DO ji = 1, npti 595 595 ! 596 596 dh_snowice(ji) = MAX( 0._wp , ( rhosn * h_s_1d(ji) + (rhoic-rau0) * h_i_1d(ji) ) * z1_rho ) … … 631 631 ! Update temperature, energy 632 632 !------------------------------------------- 633 DO ji = 1, n idx633 DO ji = 1, npti 634 634 rswitch = 1.0 - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 635 635 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 … … 637 637 638 638 DO jk = 1, nlay_s 639 DO ji = 1,n idx639 DO ji = 1,npti 640 640 ! mask enthalpy 641 641 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp, - h_s_1d(ji) ) ) … … 647 647 648 648 ! --- ensure that a_i = 0 where h_i = 0 --- 649 WHERE( h_i_1d(1:n idx) == 0._wp ) a_i_1d(1:nidx) = 0._wp649 WHERE( h_i_1d(1:npti) == 0._wp ) a_i_1d(1:npti) = 0._wp 650 650 ! 651 651 END SUBROUTINE ice_thd_dh -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90
r8564 r8565 222 222 !------------------------------------- 223 223 ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice 224 n idx = 0 ; idxice(:) = 0224 npti = 0 ; nptidx(:) = 0 225 225 DO jj = 1, jpj 226 226 DO ji = 1, jpi 227 227 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN 228 n idx = nidx+ 1229 idxice( nidx) = (jj - 1) * jpi + ji228 npti = npti + 1 229 nptidx( npti ) = (jj - 1) * jpi + ji 230 230 ENDIF 231 231 END DO … … 237 237 ! If ocean gains heat do nothing. Otherwise compute new ice formation 238 238 239 IF ( n idx> 0 ) THEN240 241 CALL tab_2d_1d( n idx, idxice(1:nidx), at_i_1d(1:nidx) , at_i )242 CALL tab_3d_2d( n idx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i (:,:,:) )243 CALL tab_3d_2d( n idx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) )244 CALL tab_3d_2d( n idx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i(:,:,:) )239 IF ( npti > 0 ) THEN 240 241 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , at_i ) 242 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) ) 243 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 244 CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 245 245 DO jl = 1, jpl 246 246 DO jk = 1, nlay_i 247 CALL tab_2d_1d( n idx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )248 END DO 249 END DO 250 CALL tab_2d_1d( n idx, idxice(1:nidx), qlead_1d (1:nidx) , qlead )251 CALL tab_2d_1d( n idx, idxice(1:nidx), t_bo_1d (1:nidx) , t_bo )252 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_opw_1d(1:nidx) , sfx_opw )253 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_opw_1d(1:nidx) , wfx_opw )254 CALL tab_2d_1d( n idx, idxice(1:nidx), zh_newice (1:nidx) , ht_i_new )255 CALL tab_2d_1d( n idx, idxice(1:nidx), zvrel_1d (1:nidx) , zvrel )256 257 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_thd_1d(1:nidx) , hfx_thd )258 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_opw_1d(1:nidx) , hfx_opw )259 CALL tab_2d_1d( n idx, idxice(1:nidx), rn_amax_1d(1:nidx) , rn_amax_2d )260 CALL tab_2d_1d( n idx, idxice(1:nidx), sss_1d (1:nidx) , sss_m )247 CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 248 END DO 249 END DO 250 CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d (1:npti) , qlead ) 251 CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d (1:npti) , t_bo ) 252 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw ) 253 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw ) 254 CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new ) 255 CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d (1:npti) , zvrel ) 256 257 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd ) 258 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw ) 259 CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d ) 260 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d (1:npti) , sss_m ) 261 261 262 262 !------------------------------------------------------------------------------| … … 265 265 DO jl = 1, jpl 266 266 DO jk = 1, nlay_i 267 WHERE( v_i_2d(1:n idx,jl) > 0._wp )268 ze_i_2d(1:n idx,jk,jl) = ze_i_2d(1:nidx,jk,jl) / v_i_2d(1:nidx,jl) * REAL( nlay_i )267 WHERE( v_i_2d(1:npti,jl) > 0._wp ) 268 ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) / v_i_2d(1:npti,jl) * REAL( nlay_i ) 269 269 ELSEWHERE 270 ze_i_2d(1:n idx,jk,jl) = 0._wp270 ze_i_2d(1:npti,jk,jl) = 0._wp 271 271 END WHERE 272 272 END DO … … 279 279 ! Keep old ice areas and volume in memory 280 280 !----------------------------------------- 281 zv_b(1:n idx,:) = v_i_2d(1:nidx,:)282 za_b(1:n idx,:) = a_i_2d(1:nidx,:)281 zv_b(1:npti,:) = v_i_2d(1:npti,:) 282 za_b(1:npti,:) = a_i_2d(1:npti,:) 283 283 284 284 !---------------------- … … 287 287 SELECT CASE ( nn_icesal ) 288 288 CASE ( 1 ) ! Sice = constant 289 zs_newice(1:n idx) = rn_icesal289 zs_newice(1:npti) = rn_icesal 290 290 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 291 DO ji = 1, n idx291 DO ji = 1, npti 292 292 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) ) 293 293 END DO 294 294 CASE ( 3 ) ! Sice = F(z) [multiyear ice] 295 zs_newice(1:n idx) = 2.3295 zs_newice(1:npti) = 2.3 296 296 END SELECT 297 297 … … 300 300 !------------------------- 301 301 ! We assume that new ice is formed at the seawater freezing point 302 DO ji = 1, n idx302 DO ji = 1, npti 303 303 ztmelts = - tmut * zs_newice(ji) ! Melting point (C) 304 304 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - ( t_bo_1d(ji) - rt0 ) ) & … … 310 310 ! Age of new ice 311 311 !---------------- 312 zo_newice(1:n idx) = 0._wp312 zo_newice(1:npti) = 0._wp 313 313 314 314 !------------------- 315 315 ! Volume of new ice 316 316 !------------------- 317 DO ji = 1, n idx317 DO ji = 1, npti 318 318 319 319 zEi = - ze_newice(ji) * r1_rhoic ! specific enthalpy of forming ice [J/kg] … … 340 340 END DO 341 341 342 zv_frazb(1:n idx) = 0._wp342 zv_frazb(1:npti) = 0._wp 343 343 IF( ln_frazil ) THEN 344 344 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 345 DO ji = 1, n idx345 DO ji = 1, npti 346 346 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) ) 347 347 zfrazb = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz … … 354 354 ! Area of new ice 355 355 !----------------- 356 DO ji = 1, n idx356 DO ji = 1, npti 357 357 za_newice(ji) = zv_newice(ji) / zh_newice(ji) 358 358 END DO … … 367 367 ! If lateral ice growth gives an ice concentration gt 1, then 368 368 ! we keep the excessive volume in memory and attribute it later to bottom accretion 369 DO ji = 1, n idx369 DO ji = 1, npti 370 370 IF ( za_newice(ji) > ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN 371 371 zda_res(ji) = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) ) … … 381 381 ! find which category to fill 382 382 DO jl = 1, jpl 383 DO ji = 1, n idx383 DO ji = 1, npti 384 384 IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN 385 385 a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice(ji) … … 389 389 END DO 390 390 END DO 391 at_i_1d(1:n idx) = SUM( a_i_2d(1:nidx,:), dim=2 )391 at_i_1d(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 ) 392 392 393 393 ! Heat content 394 DO ji = 1, n idx394 DO ji = 1, npti 395 395 jl = jcat(ji) ! categroy in which new ice is put 396 396 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice … … 398 398 399 399 DO jk = 1, nlay_i 400 DO ji = 1, n idx400 DO ji = 1, npti 401 401 jl = jcat(ji) 402 402 rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) ) … … 413 413 414 414 ! for remapping 415 h_i_old (1:n idx,0:nlay_i+1) = 0._wp416 eh_i_old(1:n idx,0:nlay_i+1) = 0._wp415 h_i_old (1:npti,0:nlay_i+1) = 0._wp 416 eh_i_old(1:npti,0:nlay_i+1) = 0._wp 417 417 DO jk = 1, nlay_i 418 DO ji = 1, n idx418 DO ji = 1, npti 419 419 h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i 420 420 eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk) … … 423 423 424 424 ! new volumes including lateral/bottom accretion + residual 425 DO ji = 1, n idx425 DO ji = 1, npti 426 426 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) ) 427 427 zv_newfra = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 ) … … 433 433 ENDDO 434 434 ! --- Ice enthalpy remapping --- ! 435 CALL ice_thd_ent( ze_i_2d(1:n idx,:,jl) )435 CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) ) 436 436 ENDDO 437 437 … … 440 440 !----------------- 441 441 DO jl = 1, jpl 442 DO ji = 1, n idx442 DO ji = 1, npti 443 443 sv_i_2d(ji,jl) = sv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) ) 444 444 END DO … … 450 450 DO jl = 1, jpl 451 451 DO jk = 1, nlay_i 452 ze_i_2d(1:n idx,jk,jl) = ze_i_2d(1:nidx,jk,jl) * v_i_2d(1:nidx,jl) * r1_nlay_i452 ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) * v_i_2d(1:npti,jl) * r1_nlay_i 453 453 END DO 454 454 END DO … … 456 456 ! 7) Change 2D vectors to 1D vectors 457 457 !------------------------------------------------------------------------------! 458 CALL tab_2d_3d( n idx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i (:,:,:) )459 CALL tab_2d_3d( n idx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) )460 CALL tab_2d_3d( n idx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i(:,:,:) )458 CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) ) 459 CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 460 CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 461 461 DO jl = 1, jpl 462 462 DO jk = 1, nlay_i 463 CALL tab_1d_2d( n idx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )464 END DO 465 END DO 466 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw )467 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw )468 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd )469 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw )463 CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 464 END DO 465 END DO 466 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_opw_1d(1:npti), sfx_opw ) 467 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_opw_1d(1:npti), wfx_opw ) 468 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d(1:npti), hfx_thd ) 469 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d(1:npti), hfx_opw ) 470 470 ! 471 ENDIF ! n idx> 0471 ENDIF ! npti > 0 472 472 ! 473 473 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90
r8562 r8565 78 78 ! 1) Cumulative integral of old enthalpy * thickness and layers interfaces 79 79 !-------------------------------------------------------------------------- 80 zeh_cum0(1:n idx,0) = 0._wp81 zh_cum0 (1:n idx,0) = 0._wp80 zeh_cum0(1:npti,0) = 0._wp 81 zh_cum0 (1:npti,0) = 0._wp 82 82 DO jk0 = 1, nlay_i+2 83 DO ji = 1, n idx83 DO ji = 1, npti 84 84 zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + eh_i_old(ji,jk0-1) 85 85 zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + h_i_old (ji,jk0-1) … … 91 91 !------------------------------------ 92 92 ! new layer thickesses 93 DO ji = 1, n idx93 DO ji = 1, npti 94 94 zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) * r1_nlay_i 95 95 END DO 96 96 97 97 ! new layers interfaces 98 zh_cum1(1:n idx,0) = 0._wp98 zh_cum1(1:npti,0) = 0._wp 99 99 DO jk1 = 1, nlay_i 100 DO ji = 1, n idx100 DO ji = 1, npti 101 101 zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji) 102 102 END DO 103 103 END DO 104 104 105 zeh_cum1(1:n idx,0) = 0._wp105 zeh_cum1(1:npti,0) = 0._wp 106 106 ! new cumulative q*h => linear interpolation 107 107 DO jk0 = 1, nlay_i+1 108 108 DO jk1 = 1, nlay_i-1 109 DO ji = 1, n idx109 DO ji = 1, npti 110 110 IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN 111 111 zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1 ) ) + & … … 117 117 END DO 118 118 ! to ensure that total heat content is strictly conserved, set: 119 zeh_cum1(1:n idx,nlay_i) = zeh_cum0(1:nidx,nlay_i+2)119 zeh_cum1(1:npti,nlay_i) = zeh_cum0(1:npti,nlay_i+2) 120 120 121 121 ! new enthalpies 122 122 DO jk1 = 1, nlay_i 123 DO ji = 1, n idx123 DO ji = 1, npti 124 124 rswitch = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 125 125 qnew(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) … … 130 130 ! comment: if input h_i_old and eh_i_old are already multiplied by a_i (as in icethd_do), 131 131 ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0 132 DO ji = 1, n idx132 DO ji = 1, npti 133 133 hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice * & 134 134 & ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90
r8564 r8565 70 70 z1_time_fl = 1._wp / rn_time_fl * rdt_ice 71 71 ! 72 DO ji = 1, n idx72 DO ji = 1, npti 73 73 74 74 !--------------------------------------------------------- -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
r8564 r8565 144 144 145 145 ! --- diag error on heat diffusion - PART 1 --- ! 146 DO ji = 1, n idx146 DO ji = 1, npti 147 147 zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 148 148 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) … … 152 152 ! 1) Initialization 153 153 !------------------ 154 DO ji = 1, n idx154 DO ji = 1, npti 155 155 isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) ) ! is there snow or not 156 156 ! layer thickness … … 159 159 END DO 160 160 ! 161 WHERE( zh_i(1:n idx) >= epsi10 ) ; z1_h_i(1:nidx) = 1._wp / zh_i(1:nidx)162 ELSEWHERE ; z1_h_i(1:n idx) = 0._wp161 WHERE( zh_i(1:npti) >= epsi10 ) ; z1_h_i(1:npti) = 1._wp / zh_i(1:npti) 162 ELSEWHERE ; z1_h_i(1:npti) = 0._wp 163 163 END WHERE 164 164 165 WHERE( zh_s(1:n idx) >= epsi10 ) ; z1_h_s(1:nidx) = 1._wp / zh_s(1:nidx)166 ELSEWHERE ; z1_h_s(1:n idx) = 0._wp165 WHERE( zh_s(1:npti) >= epsi10 ) ; z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 166 ELSEWHERE ; z1_h_s(1:npti) = 0._wp 167 167 END WHERE 168 168 ! 169 169 ! temperatures 170 ztsub (1:n idx) = t_su_1d(1:nidx) ! temperature at the previous iteration171 ztsold (1:n idx,:) = t_s_1d(1:nidx,:) ! Old snow temperature172 ztiold (1:n idx,:) = t_i_1d(1:nidx,:) ! Old ice temperature173 t_su_1d(1:n idx) = MIN( t_su_1d(1:nidx), rt0 - ztsu_err ) ! necessary170 ztsub (1:npti) = t_su_1d(1:npti) ! temperature at the previous iteration 171 ztsold (1:npti,:) = t_s_1d(1:npti,:) ! Old snow temperature 172 ztiold (1:npti,:) = t_i_1d(1:npti,:) ! Old ice temperature 173 t_su_1d(1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! necessary 174 174 ! 175 175 !------------- … … 177 177 !------------- 178 178 z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 179 DO ji = 1, n idx179 DO ji = 1, npti 180 180 ! --- Computation of i0 --- ! 181 181 ! i0 describes the fraction of solar radiation which does not contribute … … 199 199 200 200 ! --- Transmission/absorption of solar radiation in the ice --- ! 201 zradtr_s(1:n idx,0) = zftrice(1:nidx)201 zradtr_s(1:npti,0) = zftrice(1:npti) 202 202 DO jk = 1, nlay_s 203 DO ji = 1, n idx203 DO ji = 1, npti 204 204 ! ! radiation transmitted below the layer-th snow layer 205 205 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) ) … … 209 209 END DO 210 210 211 zradtr_i(1:n idx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) )211 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + zftrice(1:npti) * ( 1._wp - isnow(1:npti) ) 212 212 DO jk = 1, nlay_i 213 DO ji = 1, n idx213 DO ji = 1, npti 214 214 ! ! radiation transmitted below the layer-th ice layer 215 215 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) … … 219 219 END DO 220 220 221 ftr_ice_1d(1:n idx) = zradtr_i(1:nidx,nlay_i) ! record radiation transmitted below the ice221 ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice 222 222 ! 223 223 iconv = 0 ! number of iterations … … 228 228 iconv = iconv + 1 229 229 ! 230 ztib(1:n idx,:) = t_i_1d(1:nidx,:)231 ztsb(1:n idx,:) = t_s_1d(1:nidx,:)230 ztib(1:npti,:) = t_i_1d(1:npti,:) 231 ztsb(1:npti,:) = t_s_1d(1:npti,:) 232 232 ! 233 233 !-------------------------------- … … 236 236 IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T 237 237 ! 238 DO ji = 1, n idx238 DO ji = 1, npti 239 239 ztcond_i(ji,0) = rcdic + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 240 240 ztcond_i(ji,nlay_i) = rcdic + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) 241 241 END DO 242 242 DO jk = 1, nlay_i-1 243 DO ji = 1, n idx243 DO ji = 1, npti 244 244 ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 245 245 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) … … 249 249 ELSEIF( ln_cndi_P07 ) THEN !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 250 250 ! 251 DO ji = 1, n idx251 DO ji = 1, npti 252 252 ztcond_i(ji,0) = rcdic + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 253 253 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) … … 256 256 END DO 257 257 DO jk = 1, nlay_i-1 258 DO ji = 1, n idx258 DO ji = 1, npti 259 259 ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 260 260 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) & … … 264 264 ! 265 265 ENDIF 266 ztcond_i(1:n idx,:) = MAX( zkimin, ztcond_i(1:nidx,:) )266 ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) ) 267 267 ! 268 268 !--- G(he) : enhancement of thermal conductivity in mono-category case … … 270 270 ! Used in mono-category case only to simulate an ITD implicitly 271 271 ! Fichefet and Morales Maqueda, JGR 1997 272 zghe(1:n idx) = 1._wp272 zghe(1:npti) = 1._wp 273 273 ! 274 274 SELECT CASE ( nn_monocat ) … … 277 277 278 278 zepsilon = 0.1_wp 279 DO ji = 1, n idx279 DO ji = 1, npti 280 280 zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) ! Mean sea ice thermal conductivity 281 281 zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) ! Effective thickness he (zhe) … … 292 292 !--- Snow 293 293 DO jk = 0, nlay_s-1 294 DO ji = 1, n idx294 DO ji = 1, npti 295 295 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 296 296 END DO 297 297 END DO 298 DO ji = 1, n idx! Snow-ice interface298 DO ji = 1, npti ! Snow-ice interface 299 299 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 300 300 IF( zfac > epsi10 ) THEN … … 307 307 !--- Ice 308 308 DO jk = 0, nlay_i 309 DO ji = 1, n idx309 DO ji = 1, npti 310 310 zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 311 311 END DO 312 312 END DO 313 DO ji = 1, n idx! Snow-ice interface313 DO ji = 1, npti ! Snow-ice interface 314 314 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 315 315 END DO … … 319 319 !-------------------------------------- 320 320 DO jk = 1, nlay_i 321 DO ji = 1, n idx321 DO ji = 1, npti 322 322 zcpi = cpic + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 323 323 zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) … … 326 326 327 327 DO jk = 1, nlay_s 328 DO ji = 1, n idx328 DO ji = 1, npti 329 329 zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 330 330 END DO … … 335 335 !---------------------------- 336 336 IF ( ln_dqns_i ) THEN 337 DO ji = 1, n idx337 DO ji = 1, npti 338 338 ! update of the non solar flux according to the update in T_su 339 339 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) … … 341 341 ENDIF 342 342 343 DO ji = 1, n idx343 DO ji = 1, npti 344 344 zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH) 345 345 END DO … … 354 354 355 355 !!ice interior terms (top equation has the same form as the others) 356 ztrid (1:n idx,:,:) = 0._wp357 zindterm(1:n idx,:) = 0._wp358 zindtbis(1:n idx,:) = 0._wp359 zdiagbis(1:n idx,:) = 0._wp356 ztrid (1:npti,:,:) = 0._wp 357 zindterm(1:npti,:) = 0._wp 358 zindtbis(1:npti,:) = 0._wp 359 zdiagbis(1:npti,:) = 0._wp 360 360 361 361 DO jm = nlay_s + 2, nlay_s + nlay_i 362 DO ji = 1, n idx362 DO ji = 1, npti 363 363 jk = jm - nlay_s - 1 364 364 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) … … 370 370 371 371 jm = nlay_s + nlay_i + 1 372 DO ji = 1, n idx372 DO ji = 1, npti 373 373 !!ice bottom term 374 374 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) … … 380 380 381 381 382 DO ji = 1, n idx382 DO ji = 1, npti 383 383 ! !---------------------! 384 384 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! … … 497 497 jm_maxt = 0 498 498 jm_mint = nlay_i+5 499 DO ji = 1, n idx499 DO ji = 1, npti 500 500 jm_mint = MIN(jm_min(ji),jm_mint) 501 501 jm_maxt = MAX(jm_max(ji),jm_maxt) … … 503 503 504 504 DO jk = jm_mint+1, jm_maxt 505 DO ji = 1, n idx505 DO ji = 1, npti 506 506 jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 507 507 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) … … 510 510 END DO 511 511 512 DO ji = 1, n idx512 DO ji = 1, npti 513 513 ! ice temperatures 514 514 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) … … 516 516 517 517 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 518 DO ji = 1, n idx518 DO ji = 1, npti 519 519 jk = jm - nlay_s - 1 520 520 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) … … 522 522 END DO 523 523 524 DO ji = 1, n idx524 DO ji = 1, npti 525 525 ! snow temperatures 526 526 IF( h_s_1d(ji) > 0._wp ) THEN … … 542 542 ! zdti_max is a measure of error, it has to be under zdti_bnd 543 543 zdti_max = 0._wp 544 DO ji = 1, n idx544 DO ji = 1, npti 545 545 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 546 546 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) … … 548 548 549 549 DO jk = 1, nlay_s 550 DO ji = 1, n idx550 DO ji = 1, npti 551 551 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 552 552 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) … … 555 555 556 556 DO jk = 1, nlay_i 557 DO ji = 1, n idx557 DO ji = 1, npti 558 558 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 559 559 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) … … 576 576 ! 10) Fluxes at the interfaces 577 577 !----------------------------- 578 DO ji = 1, n idx578 DO ji = 1, npti 579 579 ! ! surface ice conduction flux 580 580 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & … … 589 589 ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 590 590 IF ( ln_dqns_i ) THEN 591 DO ji = 1, n idx591 DO ji = 1, npti 592 592 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 593 593 END DO … … 597 597 ! hfx_dif = Heat flux used to warm/cool ice in W.m-2 598 598 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 599 DO ji = 1, n idx599 DO ji = 1, npti 600 600 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 601 601 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) … … 614 614 615 615 ! --- Diagnostics SIMIP --- ! 616 DO ji = 1, n idx616 DO ji = 1, npti 617 617 !--- Conduction fluxes (positive downwards) 618 618 diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) … … 645 645 ! 646 646 DO jk = 1, nlay_i ! Sea ice energy of melting 647 DO ji = 1, n idx647 DO ji = 1, npti 648 648 ztmelts = - tmut * sz_i_1d(ji,jk) 649 649 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point … … 655 655 END DO 656 656 DO jk = 1, nlay_s ! Snow energy of melting 657 DO ji = 1, n idx657 DO ji = 1, npti 658 658 e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 659 659 END DO -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
r8564 r8565 378 378 CASE( 1 ) ! constant salinity in time and space ! 379 379 ! !---------------------------------------! 380 sz_i_1d(1:n idx,:) = rn_icesal380 sz_i_1d(1:npti,:) = rn_icesal 381 381 ! 382 382 ! !---------------------------------------------! … … 387 387 ! 388 388 ! ! Slope of the linear profile 389 WHERE( h_i_1d(1:n idx) > epsi20 ) ; z_slope_s(1:nidx) = 2._wp * s_i_1d(1:nidx) / h_i_1d(1:nidx)390 ELSEWHERE ; z_slope_s(1:n idx) = 0._wp389 WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti) 390 ELSEWHERE ; z_slope_s(1:npti) = 0._wp 391 391 END WHERE 392 392 393 393 z1_dS = 1._wp / ( zsi1 - zsi0 ) 394 DO ji = 1, n idx394 DO ji = 1, npti 395 395 zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) ) 396 396 ! ! force a constant profile when SSS too low (Baltic Sea) … … 400 400 ! Computation of the profile 401 401 DO jk = 1, nlay_i 402 DO ji = 1, n idx402 DO ji = 1, npti 403 403 ! ! linear profile with 0 surface value 404 404 zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i … … 414 414 ! !-------------------------------------------! (mean = 2.30) 415 415 ! 416 s_i_1d(1:n idx) = 2.30_wp416 s_i_1d(1:npti) = 2.30_wp 417 417 ! 418 418 !!gm cf remark in ice_var_salprof routine, CASE( 3 ) … … 420 420 zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 421 421 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) ) 422 DO ji = 1, n idx422 DO ji = 1, npti 423 423 sz_i_1d(ji,jk) = zsal 424 424 END DO
Note: See TracChangeset
for help on using the changeset viewer.
