Changeset 1572 for trunk/NEMO/LIM_SRC_3
- Timestamp:
- 2009-08-03T16:53:15+02:00 (15 years ago)
- Location:
- trunk/NEMO/LIM_SRC_3
- Files:
-
- 3 edited
-
limitd_me.F90 (modified) (4 diffs)
-
limthd.F90 (modified) (32 diffs)
-
limthd_dh.F90 (modified) (31 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limitd_me.F90
r1571 r1572 1 1 MODULE limitd_me 2 3 #if defined key_lim34 !!----------------------------------------------------------------------5 !! 'key_lim3' : LIM3 sea-ice model6 !!----------------------------------------------------------------------7 !!8 2 !!====================================================================== 9 3 !! *** MODULE limitd_me *** … … 11 5 !! computation of changes in g(h) 12 6 !!====================================================================== 13 !! 7 !! History : LIM ! 2006-02 (M. Vancoppenolle) Original code 8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & fsalt_rpo 14 9 !!---------------------------------------------------------------------- 15 !! * Modules used 10 #if defined key_lim3 11 !!---------------------------------------------------------------------- 12 !! 'key_lim3' : LIM3 sea-ice model 13 !!---------------------------------------------------------------------- 16 14 USE dom_ice 17 15 USE par_oce ! ocean parameters … … 150 148 !! Babko et al., JGR 2002 151 149 !! 152 !! ** History : 153 !! This routine is based on CICE code 154 !! and authors William H. Lipscomb, 155 !! and Elizabeth C. Hunke, LANL 156 !! are gratefully acknowledged 157 !! 158 !! (02-2006) Martin Vancoppenolle, UCL-ASTR 159 !! 150 !! This routine is based on CICE code and authors William H. Lipscomb, 151 !! and Elizabeth C. Hunke, LANL are gratefully acknowledged 160 152 !!--------------------------------------------------------------------! 161 153 !! * Arguments … … 1366 1358 !------------- 1367 1359 smsw(ji,jj) = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0 ! salt content of water frozen in voids 1368 1360 1369 1361 zsrdg2 = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge 1370 1362 1371 1363 srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity 1372 1364 1373 1365 ! ! excess of salt is flushed into the ocean 1374 fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic / rdt_ice 1375 1366 fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic / rdt_ice 1367 1376 1368 !------------------------------------ 1377 1369 ! 3.6 Increment ridging diagnostics -
trunk/NEMO/LIM_SRC_3/limthd.F90
r1571 r1572 2 2 !!====================================================================== 3 3 !! *** MODULE limthd *** 4 !! LIM thermo ice model :ice thermodynamic4 !! LIM-3 : ice thermodynamic 5 5 !!====================================================================== 6 !! History : LIM ! 2000-01 (M.A. Morales Maqueda, H. Goosse, T. Fichefet) LIM-1 7 !! 2.0 ! 2002-07 (C. Ethe, G. Madec) LIM-2 (F90 rewriting) 8 !! 3.0 ! 2005-11 (M. Vancoppenolle) LIM-3 : Multi-layer thermodynamics + salinity variations 9 !! - ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec and lim_thd_con_dif 10 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif 11 !!---------------------------------------------------------------------- 6 12 #if defined key_lim3 7 13 !!---------------------------------------------------------------------- … … 11 17 !! lim_thd_init : initialisation of sea-ice thermodynamic 12 18 !!---------------------------------------------------------------------- 13 !! * Modules used14 19 USE phycst ! physical constants 15 20 USE dom_oce ! ocean space and time domain variables 16 USE domvvl ! Variable volume17 USE lbclnk18 USE in_out_manager ! I/O manager19 21 USE ice ! LIM sea-ice variables 22 USE par_ice 20 23 USE sbc_oce ! Surface boundary condition: ocean fields 21 24 USE sbc_ice ! Surface boundary condition: ice fields 22 25 USE thd_ice ! LIM thermodynamic sea-ice variables 23 26 USE dom_ice ! LIM sea-ice domain 27 USE domvvl ! Variable volume 24 28 USE iceini 25 29 USE limthd_dif … … 28 32 USE limthd_ent 29 33 USE limtab 30 USE par_ice31 34 USE limvar 35 USE in_out_manager ! I/O manager 32 36 USE prtctl ! Print control 37 USE lbclnk 33 38 USE lib_mpp 34 39 … … 36 41 PRIVATE 37 42 38 !! * Routine accessibility 39 PUBLIC lim_thd ! called by lim_step 40 41 !! * Module variables 42 REAL(wp) :: & ! constant values 43 epsi20 = 1e-20 , & 44 epsi16 = 1e-16 , & 45 epsi06 = 1e-06 , & 46 epsi04 = 1e-04 , & 47 zzero = 0.e0 , & 48 zone = 1.e0 43 PUBLIC lim_thd ! called by lim_step 44 45 REAL(wp) :: epsi20 = 1e-20 ! constant values 46 REAL(wp) :: epsi16 = 1e-16 ! 47 REAL(wp) :: epsi06 = 1e-06 ! 48 REAL(wp) :: epsi04 = 1e-04 ! 49 REAL(wp) :: zzero = 0.e0 ! 50 REAL(wp) :: zone = 1.e0 ! 49 51 50 52 !! * Substitutions … … 52 54 # include "vectopt_loop_substitute.h90" 53 55 !!---------------------------------------------------------------------- 54 !! LIM 3.0, UCL-LOCEAN-IPSL (2005)56 !! NEMO/LIM 3.2 , UCL-LOCEAN-IPSL (2009) 55 57 !! $Id$ 56 58 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 75 77 !! - back to the geographic grid 76 78 !! 77 !! ** References : 78 !! H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90 79 !! ** References : H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90 80 !!--------------------------------------------------------------------- 81 INTEGER, INTENT(in) :: kt ! number of iteration 79 82 !! 80 !! History : 81 !! 1.0 ! 00-01 (M.A. Morales Maqueda, H. Goosse, T. Fichefet) 82 !! 2.0 ! 02-07 (C. Ethe, G. Madec) F90 83 !! 3.0 ! 05-11 (M. Vancoppenolle ) Multi-layer thermodynamics, 84 !! salinity variations 85 !!--------------------------------------------------------------------- 86 INTEGER, INTENT(in) :: kt ! number of iteration 87 !! * Local variables 88 INTEGER :: ji, jj, jk, jl, nbpb ! nb of icy pts for thermo. cal. 89 90 REAL(wp) :: & 91 zfric_umin = 5e-03 , & ! lower bound for the friction velocity 92 zfric_umax = 2e-02 ! upper bound for the friction velocity 93 94 REAL(wp) :: & 95 zinda , & ! switch for test. the val. of concen. 96 zindb, & ! switches for test. the val of arg 97 zthsnice , & 98 zfric_u , & ! friction velocity 99 zfnsol , & ! total non solar heat 100 zfontn , & ! heat flux from snow thickness 101 zfntlat, zpareff , & ! test. the val. of lead heat budget 102 zeps 103 104 REAL(wp) :: & 105 zareamin 106 83 INTEGER :: ji, jj, jk, jl ! dummy loop indices 84 INTEGER :: nbpb ! nb of icy pts for thermo. cal. 85 REAL(wp) :: zfric_umin = 5e-03 ! lower bound for the friction velocity 86 REAL(wp) :: zfric_umax = 2e-02 ! upper bound for the friction velocity 87 REAL(wp) :: zinda, zindb, zthsnice, zfric_u ! temporary scalar 88 REAL(wp) :: zfnsol, zfontn, zfntlat, zpareff ! - - 89 REAL(wp) :: zeps, zareamin, zcoef 107 90 REAL(wp), DIMENSION(jpi,jpj) :: zqlbsbq ! link with lead energy budget qldif 108 109 91 !!------------------------------------------------------------------- 110 92 111 IF( numit == nstart ) CALL lim_thd_init ! Initialization (first time-step only) 112 113 IF( kt == nit000 .AND. lwp ) THEN 114 WRITE(numout,*) 'limthd : Ice Thermodynamics' 115 WRITE(numout,*) '~~~~~~' 116 ENDIF 117 118 IF( numit == nstart ) CALL lim_thd_sal_init ! Initialization (first time-step only) 93 IF( numit == nstart ) CALL lim_thd_init ! Initialization (first time-step only) 94 95 IF( numit == nstart ) CALL lim_thd_sal_init ! Initialization (first time-step only) 96 119 97 !------------------------------------------------------------------------------! 120 98 ! 1) Initialization of diagnostic variables ! … … 125 103 ! 1.2) Heat content 126 104 !-------------------- 127 ! Change the units of heat content; from global units to 128 ! J.m3 105 ! Change the units of heat content; from global units to J.m3 129 106 130 107 DO jl = 1, jpl … … 133 110 DO ji = 1, jpi 134 111 !Energy of melting q(S,T) [J.m-3] 135 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / area(ji,jj) / & 136 MAX( v_i(ji,jj,jl) , epsi06 ) * nlay_i 112 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * nlay_i 137 113 !0 if no ice and 1 if yes 138 zindb 114 zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i(ji,jj,jl) ) ) 139 115 !convert units ! very important that this line is here 140 116 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb … … 142 118 END DO 143 119 END DO 144 END DO145 146 DO jl = 1, jpl147 120 DO jk = 1, nlay_s 148 121 DO jj = 1, jpj 149 122 DO ji = 1, jpi 150 123 !Energy of melting q(S,T) [J.m-3] 151 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / area(ji,jj) / & 152 MAX( v_s(ji,jj,jl) , epsi06 ) * nlay_s 124 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * nlay_s 153 125 !0 if no ice and 1 if yes 154 zindb 126 zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s(ji,jj,jl) ) ) 155 127 !convert units ! very important that this line is here 156 128 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb … … 180 152 ! 1.4) Compute global heat content 181 153 !----------------------------------- 182 qt_i_in (:,:)= 0.e0183 qt_s_in (:,:)= 0.e0184 qt_i_fin (:,:)= 0.e0185 qt_s_fin (:,:)= 0.e0154 qt_i_in (:,:) = 0.e0 155 qt_s_in (:,:) = 0.e0 156 qt_i_fin (:,:) = 0.e0 157 qt_s_fin (:,:) = 0.e0 186 158 sum_fluxq(:,:) = 0.e0 187 fatm (:,:) = 0.e0159 fatm (:,:) = 0.e0 188 160 189 161 ! 2) Partial computation of forcing for the thermodynamic sea ice model. ! … … 220 192 zfontn = sprecip(ji,jj) * lfus ! energy of melting 221 193 zfnsol = qns(ji,jj) ! total non solar flux 222 qldif(ji,jj) = tms(ji,jj) * ( qsr(ji,jj) &223 & + zfnsol + fdtcn(ji,jj) - zfontn &224 & + ( 1.0 - zindb ) * fsbbq(ji,jj) ) &194 qldif(ji,jj) = tms(ji,jj) * ( qsr(ji,jj) & 195 & + zfnsol + fdtcn(ji,jj) - zfontn & 196 & + ( 1.0 - zindb ) * fsbbq(ji,jj) ) & 225 197 & * ( 1.0 - at_i(ji,jj) ) * rdt_ice 226 198 … … 229 201 != 1 if positive heat budget 230 202 zpareff = 1.0 - zinda * zfntlat 231 != 0 if ice and positive heat budget and 1 if one of those two is 232 !false 233 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / & 234 MAX( at_i(ji,jj) * rdt_ice , epsi16 ) 203 != 0 if ice and positive heat budget and 1 if one of those two is false 204 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( at_i(ji,jj) * rdt_ice , epsi16 ) 235 205 236 206 ! Heat budget of the lead, energy transferred from ice to ocean … … 238 208 qdtcn (ji,jj) = zpareff * qdtcn(ji,jj) 239 209 240 ! Energy needed to bring ocean surface layer until its freezing 241 ! qcmif, limflx 242 qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj,1) & 243 & * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda ) 244 245 ! calculate oceanic heat flux (limthd_dh) 210 ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 211 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda ) 212 213 ! oceanic heat flux (limthd_dh) 246 214 fbif (ji,jj) = zindb * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) 247 215 ! … … 279 247 !------------------------------------------------------------------------------! 280 248 281 IF( lk_mpp ) CALL mpp_ini_ice(nbpb)282 283 IF (nbpb > 0) THEN ! If there is no ice, do nothing.249 IF( lk_mpp ) CALL mpp_ini_ice( nbpb ) 250 251 IF( nbpb > 0 ) THEN ! If there is no ice, do nothing. 284 252 285 253 !------------------------- … … 287 255 !------------------------- 288 256 289 CALL tab_2d_1d( nbpb, at_i_b (1:nbpb) 290 CALL tab_2d_1d( nbpb, a_i_b (1:nbpb) 291 CALL tab_2d_1d( nbpb, ht_i_b (1:nbpb) 292 CALL tab_2d_1d( nbpb, ht_s_b (1:nbpb) 293 294 CALL tab_2d_1d( nbpb, t_su_b (1:nbpb) , t_su(:,:,jl), jpi, jpj, npb(1:nbpb) )295 CALL tab_2d_1d( nbpb, sm_i_b (1:nbpb) , sm_i(:,:,jl), jpi, jpj, npb(1:nbpb) )257 CALL tab_2d_1d( nbpb, at_i_b (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) ) 258 CALL tab_2d_1d( nbpb, a_i_b (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 259 CALL tab_2d_1d( nbpb, ht_i_b (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 260 CALL tab_2d_1d( nbpb, ht_s_b (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 261 262 CALL tab_2d_1d( nbpb, t_su_b (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 263 CALL tab_2d_1d( nbpb, sm_i_b (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 296 264 DO jk = 1, nlay_s 297 CALL tab_2d_1d( nbpb, t_s_b(1:nbpb,jk) , t_s(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )298 CALL tab_2d_1d( nbpb, q_s_b(1:nbpb,jk) , e_s(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )265 CALL tab_2d_1d( nbpb, t_s_b(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 266 CALL tab_2d_1d( nbpb, q_s_b(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 299 267 END DO 300 268 DO jk = 1, nlay_i 301 CALL tab_2d_1d( nbpb, t_i_b(1:nbpb,jk) , t_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )302 CALL tab_2d_1d( nbpb, q_i_b(1:nbpb,jk) , e_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )303 CALL tab_2d_1d( nbpb, s_i_b(1:nbpb,jk) , s_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )269 CALL tab_2d_1d( nbpb, t_i_b(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 270 CALL tab_2d_1d( nbpb, q_i_b(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 271 CALL tab_2d_1d( nbpb, s_i_b(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 304 272 END DO 305 273 306 CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb) 307 CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb) 308 CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb) , fr1_i0, jpi, jpj, npb(1:nbpb) )309 CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb) , fr2_i0, jpi, jpj, npb(1:nbpb) )310 CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb) 274 CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:) , jpi, jpj, npb(1:nbpb) ) 275 CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 276 CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) ) 277 CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb), fr2_i0 , jpi, jpj, npb(1:nbpb) ) 278 CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 311 279 312 280 #if ! defined key_coupled 313 CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb) 314 CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb) 281 CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 282 CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 315 283 #endif 316 284 317 CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb) 318 CALL tab_2d_1d( nbpb, t_bo_b (1:nbpb) 319 CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb) 320 CALL tab_2d_1d( nbpb, fbif_1d (1:nbpb) 321 CALL tab_2d_1d( nbpb, qldif_1d (1:nbpb) 322 CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb) 323 CALL tab_2d_1d( nbpb, rdmsnif_1d (1:nbpb) 324 CALL tab_2d_1d( nbpb, dmgwi_1d (1:nbpb) 325 CALL tab_2d_1d( nbpb, qlbbq_1d (1:nbpb) 326 327 CALL tab_2d_1d( nbpb, fseqv_1d (1:nbpb) 328 CALL tab_2d_1d( nbpb, fsbri_1d (1:nbpb) 329 CALL tab_2d_1d( nbpb, fhbri_1d (1:nbpb) 330 CALL tab_2d_1d( nbpb, fstbif_1d (1:nbpb) 331 CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb) 285 CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 286 CALL tab_2d_1d( nbpb, t_bo_b (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) ) 287 CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) ) 288 CALL tab_2d_1d( nbpb, fbif_1d (1:nbpb), fbif , jpi, jpj, npb(1:nbpb) ) 289 CALL tab_2d_1d( nbpb, qldif_1d (1:nbpb), qldif , jpi, jpj, npb(1:nbpb) ) 290 CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb), rdmicif , jpi, jpj, npb(1:nbpb) ) 291 CALL tab_2d_1d( nbpb, rdmsnif_1d (1:nbpb), rdmsnif , jpi, jpj, npb(1:nbpb) ) 292 CALL tab_2d_1d( nbpb, dmgwi_1d (1:nbpb), dmgwi , jpi, jpj, npb(1:nbpb) ) 293 CALL tab_2d_1d( nbpb, qlbbq_1d (1:nbpb), zqlbsbq , jpi, jpj, npb(1:nbpb) ) 294 295 CALL tab_2d_1d( nbpb, fseqv_1d (1:nbpb), fseqv , jpi, jpj, npb(1:nbpb) ) 296 CALL tab_2d_1d( nbpb, fsbri_1d (1:nbpb), fsbri , jpi, jpj, npb(1:nbpb) ) 297 CALL tab_2d_1d( nbpb, fhbri_1d (1:nbpb), fhbri , jpi, jpj, npb(1:nbpb) ) 298 CALL tab_2d_1d( nbpb, fstbif_1d (1:nbpb), fstric , jpi, jpj, npb(1:nbpb) ) 299 CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb), qfvbq , jpi, jpj, npb(1:nbpb) ) 332 300 333 301 !-------------------------------- … … 335 303 !-------------------------------- 336 304 337 IF ( con_i ) CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 338 IF ( con_i ) CALL lim_thd_glohec( qt_i_in , qt_s_in , & 339 q_i_layer_in , 1 , nbpb , jl ) 340 341 !---------------------------------! 342 CALL lim_thd_dif(1,nbpb,jl) ! Ice/Snow Temperature profile ! 343 !---------------------------------! 344 345 CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 346 ! compulsory for limthd_dh 347 348 IF ( con_i ) CALL lim_thd_glohec( qt_i_fin , qt_s_fin , & 349 q_i_layer_fin , 1 , nbpb , jl ) 350 IF ( con_i ) CALL lim_thd_con_dif( 1 , nbpb , jl ) 351 352 !---------------------------------! 353 CALL lim_thd_dh(1,nbpb,jl) ! Ice/Snow thickness ! 354 !---------------------------------! 355 356 !---------------------------------! 357 CALL lim_thd_ent(1,nbpb,jl) ! Ice/Snow enthalpy remapping ! 358 !---------------------------------! 359 360 !---------------------------------! 361 CALL lim_thd_sal(1,nbpb) ! Ice salinity computation ! 362 !---------------------------------! 305 IF( con_i ) CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting 306 IF( con_i ) CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl ) 307 308 ! !---------------------------------! 309 CALL lim_thd_dif( 1, nbpb, jl ) ! Ice/Snow Temperature profile ! 310 ! !---------------------------------! 311 312 CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting compulsory for limthd_dh 313 314 IF( con_i ) CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 315 IF( con_i ) CALL lim_thd_con_dif( 1 , nbpb , jl ) 316 317 ! !---------------------------------! 318 CALL lim_thd_dh( 1, nbpb, jl ) ! Ice/Snow thickness ! 319 ! !---------------------------------! 320 321 ! !---------------------------------! 322 CALL lim_thd_ent( 1, nbpb, jl ) ! Ice/Snow enthalpy remapping ! 323 ! !---------------------------------! 324 325 ! !---------------------------------! 326 CALL lim_thd_sal( 1, nbpb ) ! Ice salinity computation ! 327 ! !---------------------------------! 363 328 364 329 ! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 365 IF ( con_i ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, & 366 q_i_layer_fin , 1 , nbpb , jl ) 367 IF ( con_i ) CALL lim_thd_con_dh ( 1 , nbpb , jl ) 330 IF( con_i ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 331 IF( con_i ) CALL lim_thd_con_dh ( 1 , nbpb , jl ) 368 332 369 333 !-------------------------------- … … 371 335 !-------------------------------- 372 336 373 CALL tab_1d_2d( nbpb, at_i , npb, at_i_b 337 CALL tab_1d_2d( nbpb, at_i , npb, at_i_b(1:nbpb), jpi, jpj ) 374 338 CALL tab_1d_2d( nbpb, ht_i(:,:,jl), npb, ht_i_b(1:nbpb), jpi, jpj ) 375 339 CALL tab_1d_2d( nbpb, ht_s(:,:,jl), npb, ht_s_b(1:nbpb), jpi, jpj ) … … 389 353 END DO 390 354 391 CALL tab_1d_2d( nbpb, fstric , npb, fstbif_1d (1:nbpb), jpi, jpj )392 CALL tab_1d_2d( nbpb, qldif , npb, qldif_1d (1:nbpb), jpi, jpj )393 CALL tab_1d_2d( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb), jpi, jpj )394 CALL tab_1d_2d( nbpb, rdmicif , npb, rdmicif_1d(1:nbpb), jpi, jpj )395 CALL tab_1d_2d( nbpb, rdmsnif , npb, rdmsnif_1d(1:nbpb), jpi, jpj )396 CALL tab_1d_2d( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb), jpi, jpj )397 CALL tab_1d_2d( nbpb, rdvosif , npb, dvsbq_1d (1:nbpb), jpi, jpj )398 CALL tab_1d_2d( nbpb, rdvobif , npb, dvbbq_1d (1:nbpb), jpi, jpj )399 CALL tab_1d_2d( nbpb, fdvolif , npb, dvlbq_1d (1:nbpb), jpi, jpj )400 CALL tab_1d_2d( nbpb, rdvonif , npb, dvnbq_1d (1:nbpb), jpi, jpj )401 CALL tab_1d_2d( nbpb, fseqv , npb, fseqv_1d (1:nbpb), jpi, jpj )402 403 IF ( num_sal .EQ.2 ) THEN404 CALL tab_1d_2d( nbpb, fsbri , npb, fsbri_1d (1:nbpb), jpi, jpj )405 CALL tab_1d_2d( nbpb, fhbri , npb, fhbri_1d (1:nbpb), jpi, jpj )355 CALL tab_1d_2d( nbpb, fstric , npb, fstbif_1d (1:nbpb), jpi, jpj ) 356 CALL tab_1d_2d( nbpb, qldif , npb, qldif_1d (1:nbpb), jpi, jpj ) 357 CALL tab_1d_2d( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb), jpi, jpj ) 358 CALL tab_1d_2d( nbpb, rdmicif, npb, rdmicif_1d(1:nbpb), jpi, jpj ) 359 CALL tab_1d_2d( nbpb, rdmsnif, npb, rdmsnif_1d(1:nbpb), jpi, jpj ) 360 CALL tab_1d_2d( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb), jpi, jpj ) 361 CALL tab_1d_2d( nbpb, rdvosif, npb, dvsbq_1d (1:nbpb), jpi, jpj ) 362 CALL tab_1d_2d( nbpb, rdvobif, npb, dvbbq_1d (1:nbpb), jpi, jpj ) 363 CALL tab_1d_2d( nbpb, fdvolif, npb, dvlbq_1d (1:nbpb), jpi, jpj ) 364 CALL tab_1d_2d( nbpb, rdvonif, npb, dvnbq_1d (1:nbpb), jpi, jpj ) 365 CALL tab_1d_2d( nbpb, fseqv , npb, fseqv_1d (1:nbpb), jpi, jpj ) 366 367 IF( num_sal == 2 ) THEN 368 CALL tab_1d_2d( nbpb, fsbri, npb, fsbri_1d(1:nbpb), jpi, jpj ) 369 CALL tab_1d_2d( nbpb, fhbri, npb, fhbri_1d(1:nbpb), jpi, jpj ) 406 370 ENDIF 407 371 408 372 !+++++ 409 !temporary stuff for a dummy version373 !temporary stuff for a dummy version 410 374 CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb) , jpi, jpj ) 411 375 CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb) , jpi, jpj ) … … 417 381 !+++++ 418 382 419 IF( lk_mpp ) CALL mpp_comm_free(ncomm_ice) !RB necessary ??420 ENDIF ! nbpb421 422 END DO ! jl383 IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? 384 ENDIF 385 ! 386 END DO 423 387 424 388 !------------------------------------------------------------------------------! … … 431 395 432 396 ! Enthalpies are global variables we have to readjust the units 397 zcoef = 1.e0 / ( unit_fac * REAL(nlay_i) ) 433 398 DO jl = 1, jpl 434 399 DO jk = 1, nlay_i 435 DO jj = 1, jpj 436 DO ji = 1, jpi 437 ! Change dimensions 438 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 439 440 ! Multiply by volume, divide by nlayers so that heat content in 10^9 Joules 441 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 442 area(ji,jj) * a_i(ji,jj,jl) * & 443 ht_i(ji,jj,jl) / nlay_i 444 END DO !ji 445 END DO !jj 446 END DO !jk 447 END DO !jl 400 ! Multiply by volume, divide by nlayers so that heat content in 10^9 Joules 401 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) * zcoef 402 END DO 403 END DO 448 404 449 405 !------------------------ … … 452 408 453 409 ! Enthalpies are global variables we have to readjust the units 410 zcoef = 1.e0 / ( unit_fac * REAL(nlay_s) ) 454 411 DO jl = 1, jpl 455 412 DO jk = 1, nlay_s 456 DO jj = 1, jpj 457 DO ji = 1, jpi 458 ! Change dimensions 459 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 460 ! Multiply by volume, so that heat content in 10^9 Joules 461 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 462 a_i(ji,jj,jl) * ht_s(ji,jj,jl) / nlay_s 463 END DO !ji 464 END DO !jj 465 END DO !jk 466 END DO !jl 413 ! Multiply by volume, so that heat content in 10^9 Joules 414 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) * zcoef 415 END DO 416 END DO 467 417 468 418 !---------------------------------- … … 474 424 ! 5.4) Diagnostic thermodynamic growth rates 475 425 !-------------------------------------------- 476 d_v_i_thd (:,:,:) = v_i(:,:,:)- old_v_i(:,:,:) ! ice volumes477 dv_dt_thd(:,:,:) 478 479 IF ( con_i )fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)426 d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes 427 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) / rdt_ice * 86400.0 428 429 IF( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 480 430 481 431 IF(ln_ctl) THEN ! Control print … … 514 464 END SUBROUTINE lim_thd 515 465 516 !=============================================================================== 517 518 SUBROUTINE lim_thd_glohec(eti,ets,etilayer,kideb,kiut,jl) 466 467 SUBROUTINE lim_thd_glohec( eti, ets, etilayer, kideb, kiut, jl ) 519 468 !!----------------------------------------------------------------------- 520 469 !! *** ROUTINE lim_thd_glohec *** … … 522 471 !! ** Purpose : Compute total heat content for each category 523 472 !! Works with 1d vectors only 473 !!----------------------------------------------------------------------- 474 INTEGER , INTENT(in ) :: kideb, kiut ! bounds for the spatial loop 475 INTEGER , INTENT(in ) :: jl ! category number 476 REAL(wp), INTENT( out), DIMENSION (jpij,jpl ) :: eti, ets ! vertically-summed heat content for ice & snow 477 REAL(wp), INTENT( out), DIMENSION (jpij,jkmax) :: etilayer ! heat content for ice layers 524 478 !! 525 !! history : 526 !! 9.9 ! 07-04 (M.Vancoppenolle) original code 479 INTEGER :: ji,jk ! loop indices 480 REAL(wp) :: zdes ! snow heat content increment (dummy) 481 REAL(wp) :: zeps ! very small value (1.e-10) 527 482 !!----------------------------------------------------------------------- 528 !! * Local variables 529 INTEGER, INTENT(in) :: & 530 kideb, kiut, & ! bounds for the spatial loop 531 jl ! category number 532 533 REAL(wp), DIMENSION (jpij,jpl), INTENT(out) :: & 534 eti, ets ! vertically-summed heat content for ice /snow 535 536 REAL(wp), DIMENSION (jpij,jkmax), INTENT(out) :: & 537 etilayer ! heat content for ice layers 538 539 REAL(wp) :: & 540 zdes, & ! snow heat content increment (dummy) 541 zeps ! very small value (1.e-10) 542 543 INTEGER :: & 544 ji,jk ! loop indices 545 546 !!----------------------------------------------------------------------- 547 eti(:,:) = 0.0 548 ets(:,:) = 0.0 483 eti(:,:) = 0.e0 484 ets(:,:) = 0.e0 549 485 zeps = 1.0e-10 550 486 551 ! total q over all layers, ice [J.m-2] 552 DO jk = 1, nlay_i 487 DO jk = 1, nlay_i ! total q over all layers, ice [J.m-2] 553 488 DO ji = kideb, kiut 554 etilayer(ji,jk) = q_i_b(ji,jk) & 555 * ht_i_b(ji) / nlay_i 556 eti(ji,jl) = eti(ji,jl) + etilayer(ji,jk) 489 etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 490 eti (ji,jl) = eti(ji,jl) + etilayer(ji,jk) 557 491 END DO 558 492 END DO 559 560 ! total q over all layers, snow [J.m-2] 561 DO ji = kideb, kiut 562 zdes = q_s_b(ji,1) * ht_s_b(ji) / nlay_s 563 ets(ji,jl) = ets(ji,jl) + zdes 564 END DO 565 566 WRITE(numout,*) ' lim_thd_glohec ' 567 WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) / rdt_ice 568 WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) / rdt_ice 569 WRITE(numout,*) ' qt_in : ', ( eti(jiindex_1d,jl) + & 570 ets(jiindex_1d,jl) ) / rdt_ice 571 493 DO ji = kideb, kiut ! total q over all layers, snow [J.m-2] 494 ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / nlay_s 495 END DO 496 497 IF(lwp) WRITE(numout,*) ' lim_thd_glohec ' 498 IF(lwp) WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) / rdt_ice 499 IF(lwp) WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) / rdt_ice 500 IF(lwp) WRITE(numout,*) ' qt_in : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) / rdt_ice 501 ! 572 502 END SUBROUTINE lim_thd_glohec 573 503 574 !=============================================================================== 575 576 SUBROUTINE lim_thd_con_dif(kideb,kiut,jl) 504 505 SUBROUTINE lim_thd_con_dif( kideb, kiut, jl ) 577 506 !!----------------------------------------------------------------------- 578 507 !! *** ROUTINE lim_thd_con_dif *** 579 508 !! 580 509 !! ** Purpose : Test energy conservation after heat diffusion 581 !!582 !! history :583 !! 9.9 ! 07-04 (M.Vancoppenolle) original code584 510 !!------------------------------------------------------------------- 585 !! * Local variables 586 INTEGER, INTENT(in) :: & 587 kideb, kiut, & !: bounds for the spatial loop 588 jl !: category number 589 590 REAL(wp) :: & !: ! goes to trash 591 meance, & !: mean conservation error 592 max_cons_err, & !: maximum tolerated conservation error 593 max_surf_err !: maximum tolerated surface error 594 595 INTEGER :: & 596 numce !: number of points for which conservation 597 ! is violated 598 INTEGER :: & 599 ji,jk, & !: loop indices 600 zji, zjj 511 INTEGER , INTENT(in ) :: kideb, kiut ! bounds for the spatial loop 512 INTEGER , INTENT(in ) :: jl ! category number 513 514 INTEGER :: ji, jk ! loop indices 515 INTEGER :: zji, zjj 516 INTEGER :: numce ! number of points for which conservation is violated 517 REAL(wp) :: meance ! mean conservation error 518 REAL(wp) :: max_cons_err, max_surf_err 601 519 !!--------------------------------------------------------------------- 602 520 603 max_cons_err = 1.0 604 max_surf_err = 0.001 521 max_cons_err = 1.0 ! maximum tolerated conservation error 522 max_surf_err = 0.001 ! maximum tolerated surface error 605 523 606 524 !-------------------------- … … 609 527 ! global 610 528 DO ji = kideb, kiut 611 dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) & 612 + qt_s_fin(ji,jl) - qt_s_in(ji,jl) 529 dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) + qt_s_fin(ji,jl) - qt_s_in(ji,jl) 613 530 END DO 614 531 ! layer by layer … … 620 537 621 538 DO ji = kideb, kiut 622 zji = MOD( npb(ji) - 1, jpi ) + 1 623 zjj = ( npb(ji) - 1 ) / jpi + 1 624 625 fatm(ji,jl) = & 626 qnsr_ice_1d(ji) + & ! atm non solar 627 (1.0-i0(ji))*qsr_ice_1d(ji) ! atm solar 628 629 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji)*i0(ji) & 630 - fstroc(zji,zjj,jl) 539 zji = MOD( npb(ji) - 1, jpi ) + 1 540 zjj = ( npb(ji) - 1 ) / jpi + 1 541 542 fatm(ji,jl) = qnsr_ice_1d(ji) + (1.0-i0(ji))*qsr_ice_1d(ji) 543 544 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji)*i0(ji) - fstroc(zji,zjj,jl) 631 545 END DO 632 546 … … 651 565 WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err 652 566 WRITE(numout,*) ' After lim_thd_dif, category : ', jl 653 WRITE(numout,*) ' Mean conservation error on big error points ', meance, & 654 numit 567 WRITE(numout,*) ' Mean conservation error on big error points ', meance, numit 655 568 WRITE(numout,*) ' Number of points where there is a cons err gt than c.e. : ', numce, numit 656 569 … … 751 664 752 665 END DO 753 666 ! 754 667 END SUBROUTINE lim_thd_con_dif 755 668 756 !==============================================================================757 669 758 670 SUBROUTINE lim_thd_con_dh(kideb,kiut,jl) … … 765 677 !! 9.9 ! 07-04 (M.Vancoppenolle) original code 766 678 !!----------------------------------------------------------------------- 767 !! * Local variables768 679 INTEGER, INTENT(in) :: & 769 680 kideb, kiut, & !: bounds for the spatial loop … … 882 793 883 794 ENDIF 884 885 END DO 886 795 ! 796 END DO 797 ! 887 798 END SUBROUTINE lim_thd_con_dh 888 !============================================================================== 889 890 SUBROUTINE lim_thd_enmelt( kideb,kiut)799 800 801 SUBROUTINE lim_thd_enmelt( kideb, kiut ) 891 802 !!----------------------------------------------------------------------- 892 803 !! *** ROUTINE lim_thd_enmelt *** … … 896 807 !! ** Method : Formula (Bitz and Lipscomb, 1999) 897 808 !! 898 !! history : Martin Vancoppenolle, May 2007899 809 !!------------------------------------------------------------------- 900 INTEGER, INTENT(in) :: & 901 kideb, kiut !: bounds for the spatial loop 902 903 REAL(wp) :: & !: goes to trash 904 ztmelts , & !: sea ice freezing point in K 905 zeps 906 907 INTEGER :: & 908 ji, & !: spatial loop index 909 jk !: vertical index 910 810 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 811 !! 812 INTEGER :: ji, jk !dummy loop indices 813 REAL(wp) :: ztmelts, zeps ! temporary scalar 911 814 !!------------------------------------------------------------------- 912 zeps = 1.0e-10 913 914 ! Sea ice energy of melting 915 DO jk = 1, nlay_i 815 zeps = 1.e-10 816 ! 817 DO jk = 1, nlay_i ! Sea ice energy of melting 916 818 DO ji = kideb, kiut 917 ztmelts = - tmut * s_i_b(ji,jk) + rtt 918 q_i_b(ji,jk) = rhoic*( cpic * ( ztmelts - t_i_b(ji,jk) ) & 919 + lfus * ( 1.0 - (ztmelts-rtt)/MIN((t_i_b(ji,jk)-rtt),-zeps) ) & 920 - rcp * ( ztmelts-rtt ) ) 921 END DO !ji 922 END DO !jk 923 924 ! Snow energy of melting 925 DO jk = 1, nlay_s 819 ztmelts = - tmut * s_i_b(ji,jk) + rtt 820 q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) & 821 & + lfus * ( 1.0 - (ztmelts-rtt) / MIN( t_i_b(ji,jk)-rtt, -zeps ) ) & 822 & - rcp * ( ztmelts-rtt ) ) 823 END DO 824 END DO 825 DO jk = 1, nlay_s ! Snow energy of melting 926 826 DO ji = kideb,kiut 927 827 q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 928 END DO !ji929 END DO !jk930 828 END DO 829 END DO 830 ! 931 831 END SUBROUTINE lim_thd_enmelt 932 832 933 !==============================================================================934 833 935 834 SUBROUTINE lim_thd_init … … 939 838 !! 940 839 !! ** Purpose : Physical constants and parameters linked to the ice 941 !! thermodynamics840 !! thermodynamics 942 841 !! 943 842 !! ** Method : Read the namicethd namelist and check the ice-thermo 944 !! parameter values called at the first timestep (nit000)843 !! parameter values called at the first timestep (nit000) 945 844 !! 946 845 !! ** input : Namelist namicether 947 !! 948 !! history : 949 !! 8.5 ! 03-08 (C. Ethe) original code 950 !!------------------------------------------------------------------- 951 NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, & 952 & hicmin, hiclim, amax , & 953 & sbeta , parlat, hakspl, hibspl, exld, & 954 & hakdif, hnzst , thth , parsub, alphs, betas, & 846 !!------------------------------------------------------------------- 847 NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, & 848 & hicmin, hiclim, amax , & 849 & sbeta , parlat, hakspl, hibspl, exld, & 850 & hakdif, hnzst , thth , parsub, alphs, betas, & 955 851 & kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi 956 852 !!------------------------------------------------------------------- 957 853 958 ! Define the initial parameters 959 ! ------------------------- 960 REWIND( numnam_ice ) 854 IF(lwp) THEN 855 WRITE(numout,*) 856 WRITE(numout,*) 'lim_thd : Ice Thermodynamics' 857 WRITE(numout,*) '~~~~~~~' 858 ENDIF 859 860 REWIND( numnam_ice ) ! read Namelist numnam_ice 961 861 READ ( numnam_ice , namicethd ) 962 IF (lwp) THEN 862 863 IF(lwp) THEN ! control print 963 864 WRITE(numout,*) 964 WRITE(numout,*)'lim_thd_init : ice parameters for ice thermodynamic computation ' 965 WRITE(numout,*)'~~~~~~~~~~~~' 966 WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt 967 WRITE(numout,*)' ice thick. for lateral accretion in NH (SH) hiccrit(1/2) = ', hiccrit 968 WRITE(numout,*)' Frazil ice thickness as a function of wind or not fraz_swi = ', fraz_swi 969 WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom maxfrazb = ', maxfrazb 970 WRITE(numout,*)' Thresold relative drift speed for collection of frazil vfrazb = ', vfrazb 971 WRITE(numout,*)' Squeezing coefficient for collection of frazil Cfrazb = ', Cfrazb 972 WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin 973 WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim 974 WRITE(numout,*)' maximum lead fraction amax = ', amax 975 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 976 WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta 977 WRITE(numout,*)' percentage of energy used for lateral ablation parlat = ', parlat 978 WRITE(numout,*)' slope of distr. for Hakkinen-Mellor lateral melting hakspl = ', hakspl 979 WRITE(numout,*)' slope of distribution for Hibler lateral melting hibspl = ', hibspl 980 WRITE(numout,*)' exponent for leads-closure rate exld = ', exld 981 WRITE(numout,*)' coefficient for diffusions of ice and snow hakdif = ', hakdif 982 WRITE(numout,*)' threshold thick. for comp. of eq. thermal conductivity zhth = ', thth 983 WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst 984 WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub 985 WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs 986 WRITE(numout,*)' coefficient for ice-lead partition of snowfall betas = ', betas 987 WRITE(numout,*)' extinction radiation parameter in sea ice (1.0) kappa_i = ', kappa_i 988 WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nconv_i_thd = ', nconv_i_thd 989 WRITE(numout,*)' maximal err. on T for heat diffusion computation maxer_i_thd = ', maxer_i_thd 990 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice thcon_i_swi = ', thcon_i_swi 991 WRITE(numout,*) 865 WRITE(numout,*)' Namelist of ice parameters for ice thermodynamic computation ' 866 WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt 867 WRITE(numout,*)' ice thick. for lateral accretion in NH (SH) hiccrit(1/2) = ', hiccrit 868 WRITE(numout,*)' Frazil ice thickness as a function of wind or not fraz_swi = ', fraz_swi 869 WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom maxfrazb = ', maxfrazb 870 WRITE(numout,*)' Thresold relative drift speed for collection of frazil vfrazb = ', vfrazb 871 WRITE(numout,*)' Squeezing coefficient for collection of frazil Cfrazb = ', Cfrazb 872 WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin 873 WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim 874 WRITE(numout,*)' maximum lead fraction amax = ', amax 875 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 876 WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta 877 WRITE(numout,*)' percentage of energy used for lateral ablation parlat = ', parlat 878 WRITE(numout,*)' slope of distr. for Hakkinen-Mellor lateral melting hakspl = ', hakspl 879 WRITE(numout,*)' slope of distribution for Hibler lateral melting hibspl = ', hibspl 880 WRITE(numout,*)' exponent for leads-closure rate exld = ', exld 881 WRITE(numout,*)' coefficient for diffusions of ice and snow hakdif = ', hakdif 882 WRITE(numout,*)' threshold thick. for comp. of eq. thermal conductivity zhth = ', thth 883 WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst 884 WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub 885 WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs 886 WRITE(numout,*)' coefficient for ice-lead partition of snowfall betas = ', betas 887 WRITE(numout,*)' extinction radiation parameter in sea ice (1.0) kappa_i = ', kappa_i 888 WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nconv_i_thd = ', nconv_i_thd 889 WRITE(numout,*)' maximal err. on T for heat diffusion computation maxer_i_thd = ', maxer_i_thd 890 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice thcon_i_swi = ', thcon_i_swi 992 891 ENDIF 993 892 ! 994 893 rcdsn = hakdif * rcdsn 995 894 rcdic = hakdif * rcdic 996 997 895 ! 998 896 END SUBROUTINE lim_thd_init 999 897 1000 898 #else 1001 !!====================================================================== 1002 !! *** MODULE limthd *** 1003 !! No sea ice model 1004 !!====================================================================== 899 !!---------------------------------------------------------------------- 900 !! Default option NO LIM3 sea-ice model 901 !!---------------------------------------------------------------------- 1005 902 CONTAINS 1006 903 SUBROUTINE lim_thd ! Empty routine -
trunk/NEMO/LIM_SRC_3/limthd_dh.F90
r1571 r1572 1 1 MODULE limthd_dh 2 !!====================================================================== 3 !! *** MODULE limthd_dh *** 4 !! LIM-3 : thermodynamic growth and decay of the ice 5 !!====================================================================== 6 !! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D 7 !! ! 2005-06 (M. Vancoppenolle) 3D version 8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif & rdmicif 9 !!---------------------------------------------------------------------- 2 10 #if defined key_lim3 3 11 !!---------------------------------------------------------------------- 4 12 !! 'key_lim3' LIM3 sea-ice model 5 13 !!---------------------------------------------------------------------- 6 !!====================================================================== 7 !! *** MODULE limthd_dh *** 8 !! thermodynamic growth and decay of the ice 9 !!====================================================================== 10 14 !! lim_thd_dh : vertical accr./abl. and lateral ablation of sea ice 11 15 !!---------------------------------------------------------------------- 12 !! lim_thd_dh : vertical accr./abl. and lateral ablation of sea ice13 !! * Modules used14 15 16 USE par_oce ! ocean parameters 16 17 USE phycst ! physical constants (OCE directory) … … 27 28 PRIVATE 28 29 29 !! * Routine accessibility 30 PUBLIC lim_thd_dh ! called by lim_thd 31 32 !! * Module variables 33 REAL(wp) :: & ! constant values 34 epsi20 = 1e-20 , & 35 epsi13 = 1e-13 , & 36 epsi16 = 1e-16 , & 37 zzero = 0.e0 , & 38 zone = 1.e0 30 PUBLIC lim_thd_dh ! called by lim_thd 31 32 REAL(wp) :: epsi20 = 1e-20 ! constant values 33 REAL(wp) :: epsi13 = 1e-13 ! 34 REAL(wp) :: epsi16 = 1e-16 ! 35 REAL(wp) :: zzero = 0.e0 ! 36 REAL(wp) :: zone = 1.e0 ! 39 37 40 38 !!---------------------------------------------------------------------- 41 !! LIM 3.0, UCL-LOCEAN-IPSL (2005)39 !! NEMO/LIM 3.2, UCL-LOCEAN-IPSL (2009) 42 40 !! $Id$ 43 41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 49 47 !!------------------------------------------------------------------ 50 48 !! *** ROUTINE lim_thd_dh *** 49 !! 50 !! ** Purpose : determines variations of ice and snow thicknesses. 51 !! 52 !! ** Method : Ice/Snow surface melting arises from imbalance in surface fluxes 53 !! Bottom accretion/ablation arises from flux budget 54 !! Snow thickness can increase by precipitation and decrease by sublimation 55 !! If snow load excesses Archmiede limit, snow-ice is formed by 56 !! the flooding of sea-water in the snow 57 !! 58 !! 1) Compute available flux of heat for surface ablation 59 !! 2) Compute snow and sea ice enthalpies 60 !! 3) Surface ablation and sublimation 61 !! 4) Bottom accretion/ablation 62 !! 5) Case of Total ablation 63 !! 6) Snow ice formation 64 !! 65 !! References : Bitz and Lipscomb, 1999, J. Geophys. Res. 66 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 67 !! Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let. 68 !! Vancoppenolle et al.,2009, Ocean Modelling 51 69 !!------------------------------------------------------------------ 52 !! ** Purpose : 53 !! This routine determines variations of ice and snow thicknesses. 54 !! ** Method : 55 !! Ice/Snow surface melting arises from imbalance in surface fluxes 56 !! Bottom accretion/ablation arises from flux budget 57 !! Snow thickness can increase by precipitation and decrease by 58 !! sublimation 59 !! If snow load excesses Archmiede limit, snow-ice is formed by 60 !! the flooding of sea-water in the snow 61 !! ** Steps 62 !! 1) Compute available flux of heat for surface ablation 63 !! 2) Compute snow and sea ice enthalpies 64 !! 3) Surface ablation and sublimation 65 !! 4) Bottom accretion/ablation 66 !! 5) Case of Total ablation 67 !! 6) Snow ice formation 68 !! 69 !! ** Arguments 70 !! 71 !! ** Inputs / Outputs 72 !! 73 !! ** External 74 !! 75 !! ** References : Bitz and Lipscomb, JGR 99 76 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 77 !! Vancoppenolle, Fichefet and Bitz, GRL 2005 78 !! Vancoppenolle et al., OM08 79 !! 80 !! ** History : 81 !! original code 01-04 (LIM) 82 !! original routine 83 !! (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium 84 !! (06/07-2005) 3D version 85 !! (03-2008) Clean code 86 !! 70 INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied 71 INTEGER , INTENT(in) :: jl ! Thickness cateogry number 72 !! 73 INTEGER :: ji , jk ! dummy loop indices 74 INTEGER :: zji, zjj ! 2D corresponding indices to ji 75 INTEGER :: isnow ! switch for presence (1) or absence (0) of snow 76 INTEGER :: isnowic ! snow ice formation not 77 INTEGER :: i_ice_switch ! ice thickness above a certain treshold or not 78 INTEGER :: iter 79 80 REAL(wp) :: zzfmass_i, zzfmass_s ! temporary scalar 81 REAL(wp) :: zhsnew, zihgnew, ztmelts ! temporary scalar 82 REAL(wp) :: zhn, zdhcf, zdhbf, zhni, zhnfi, zihg ! 83 REAL(wp) :: zdhnm, zhnnew, zhisn, zihic ! 84 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 85 REAL(wp) :: zds ! increment of bottom ice salinity 86 REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads 87 REAL(wp) :: zsm_snowice ! snow-ice salinity 88 REAL(wp) :: zswi1 ! switch for computation of bottom salinity 89 REAL(wp) :: zswi12 ! switch for computation of bottom salinity 90 REAL(wp) :: zswi2 ! switch for computation of bottom salinity 91 REAL(wp) :: zgrr ! bottom growth rate 92 REAL(wp) :: ztform ! bottom formation temperature 93 94 REAL(wp), DIMENSION(jpij) :: zh_i ! ice layer thickness 95 REAL(wp), DIMENSION(jpij) :: zh_s ! snow layer thickness 96 REAL(wp), DIMENSION(jpij) :: ztfs ! melting point 97 REAL(wp), DIMENSION(jpij) :: zhsold ! old snow thickness 98 REAL(wp), DIMENSION(jpij) :: zqprec ! energy of fallen snow 99 REAL(wp), DIMENSION(jpij) :: zqfont_su ! incoming, remaining surface energy 100 REAL(wp), DIMENSION(jpij) :: zqfont_bo ! incoming, bottom energy 101 REAL(wp), DIMENSION(jpij) :: z_f_surf ! surface heat for ablation 102 REAL(wp), DIMENSION(jpij) :: zhgnew ! new ice thickness 103 REAL(wp), DIMENSION(jpij) :: zfmass_i ! 104 105 REAL(wp), DIMENSION(jpij) :: zdh_s_mel ! snow melt 106 REAL(wp), DIMENSION(jpij) :: zdh_s_pre ! snow precipitation 107 REAL(wp), DIMENSION(jpij) :: zdh_s_sub ! snow sublimation 108 REAL(wp), DIMENSION(jpij) :: zfsalt_melt ! salt flux due to ice melt 109 110 REAL(wp) , DIMENSION(jpij,jkmax) :: zdeltah 111 112 ! Pathological cases 113 REAL(wp), DIMENSION(jpij) :: zfdt_init ! total incoming heat for ice melt 114 REAL(wp), DIMENSION(jpij) :: zfdt_final ! total remaing heat for ice melt 115 REAL(wp), DIMENSION(jpij) :: zqt_i ! total ice heat content 116 REAL(wp), DIMENSION(jpij) :: zqt_s ! total snow heat content 117 REAL(wp), DIMENSION(jpij) :: zqt_dummy ! dummy heat content 118 119 REAL(wp), DIMENSION(jpij,jkmax) :: zqt_i_lay ! total ice heat content 120 121 ! Heat conservation 122 INTEGER :: num_iter_max, numce_dh 123 REAL(wp) :: meance_dh 124 INTEGER , DIMENSION(jpij) :: innermelt 125 REAL(wp), DIMENSION(jpij) :: zfbase, zdq_i 87 126 !!------------------------------------------------------------------ 88 !! * Arguments89 INTEGER , INTENT (IN) :: &90 kideb , & !: Start point on which the the computation is applied91 kiut , & !: End point on which the the computation is applied92 jl !: Thickness cateogry number93 94 !! * Local variables95 INTEGER :: &96 ji , & !: space index97 jk , & !: ice layer index98 isnow , & !: switch for presence (1) or absence (0) of snow99 zji , & !: 2D corresponding indices to ji100 zjj , &101 isnowic , & !: snow ice formation not102 i_ice_switch , & !: ice thickness above a certain treshold or not103 iter104 105 REAL(wp) :: &106 zzfmass_i , &107 zzfmass_s , &108 zhsnew , & !: new snow thickness109 zihgnew , & !: switch for total ablation110 ztmelts , & !: melting point111 zhn , &112 zdhcf , &113 zdhbf , &114 zhni , &115 zhnfi , &116 zihg , &117 zdhnm , &118 zhnnew , &119 zeps = 1.0e-13, &120 zhisn , &121 zfracs , & !: fractionation coefficient for bottom salt122 !: entrapment123 zds , & !: increment of bottom ice salinity124 zcoeff , & !: dummy argument for snowfall partitioning125 !: over ice and leads126 zsm_snowice, & !: snow-ice salinity127 zswi1 , & !: switch for computation of bottom salinity128 zswi12 , & !: switch for computation of bottom salinity129 zswi2 , & !: switch for computation of bottom salinity130 zgrr , & !: bottom growth rate131 zihic , & !:132 ztform !: bottom formation temperature133 134 REAL(wp) , DIMENSION(jpij) :: &135 zh_i , & ! ice layer thickness136 zh_s , & ! snow layer thickness137 ztfs , & ! melting point138 zhsold , & ! old snow thickness139 zqprec , & !: energy of fallen snow140 zqfont_su , & ! incoming, remaining surface energy141 zqfont_bo ! incoming, bottom energy142 143 REAL(wp) , DIMENSION(jpij) :: &144 z_f_surf, & ! surface heat for ablation145 zhgnew , & ! new ice thickness146 zfmass_i147 148 REAL(wp), DIMENSION(jpij) :: &149 zdh_s_mel , & ! snow melt150 zdh_s_pre , & ! snow precipitation151 zdh_s_sub , & ! snow sublimation152 zfsalt_melt ! salt flux due to ice melt153 154 REAL(wp) , DIMENSION(jpij,jkmax) :: &155 zdeltah156 157 ! Pathological cases158 REAL(wp), DIMENSION(jpij) :: &159 zfdt_init , & !: total incoming heat for ice melt160 zfdt_final , & !: total remaing heat for ice melt161 zqt_i , & !: total ice heat content162 zqt_s , & !: total snow heat content163 zqt_dummy !: dummy heat content164 165 REAL(wp), DIMENSION(jpij,jkmax) :: &166 zqt_i_lay !: total ice heat content167 168 ! Heat conservation169 REAL(wp), DIMENSION(jpij) :: &170 zfbase, &171 zdq_i172 173 INTEGER, DIMENSION(jpij) :: &174 innermelt175 176 REAL(wp) :: &177 meance_dh178 179 INTEGER :: &180 num_iter_max, &181 numce_dh182 127 183 128 zfsalt_melt(:) = 0.0 … … 198 143 isnow = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) ) 199 144 ztfs(ji) = isnow * rtt + ( 1.0 - isnow ) * rtt 200 z_f_surf(ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * & 201 qsr_ice_1d(ji) - fc_su(ji) 202 z_f_surf(ji) = MAX( zzero , z_f_surf(ji) ) * & 203 MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) ) 204 zfdt_init(ji) = ( z_f_surf(ji) + & 205 MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) & 206 * rdt_ice 145 z_f_surf(ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 146 z_f_surf(ji) = MAX( zzero , z_f_surf(ji) ) * MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) ) 147 zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice 207 148 END DO ! ji 208 149 … … 226 167 DO jk = 1, nlay_s 227 168 DO ji = kideb,kiut 228 zqt_s(ji) 169 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s 229 170 END DO 230 171 END DO … … 235 176 DO ji = kideb,kiut 236 177 zqt_i(ji) = zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 237 zqt_i_lay(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i178 zqt_i_lay(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 238 179 END DO 239 180 END DO … … 267 208 DO ji = kideb, kiut 268 209 ! tatm_ice is now in K 269 zqprec(ji) = rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) 270 zqfont_su(ji) = z_f_surf(ji) * rdt_ice 271 zdeltah(ji,1) = MIN( 0.0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) ) 272 zqfont_su(ji) = MAX( 0.0 , - zdh_s_pre(ji) - zdeltah(ji,1) ) * & 273 zqprec(ji) 274 zdeltah(ji,1) = MAX( - zdh_s_pre(ji) , zdeltah(ji,1) ) 275 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) 210 zqprec (ji) = rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) 211 zqfont_su(ji) = z_f_surf(ji) * rdt_ice 212 zdeltah (ji,1) = MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) ) 213 zqfont_su(ji) = MAX( 0.e0 , - zdh_s_pre(ji) - zdeltah(ji,1) ) * zqprec(ji) 214 zdeltah (ji,1) = MAX( - zdh_s_pre(ji) , zdeltah(ji,1) ) 215 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) 276 216 ! heat conservation 277 qt_s_in(ji,jl) = qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji)278 zqt_s (ji) = zqt_s(ji)+ zqprec(ji) * zdh_s_pre(ji)279 zqt_s (ji) = MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 )217 qt_s_in(ji,jl) = qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji) 218 zqt_s (ji) = zqt_s (ji) + zqprec(ji) * zdh_s_pre(ji) 219 zqt_s (ji) = MAX( zqt_s(ji) - zqfont_su(ji) , 0.e0 ) 280 220 END DO 281 221 … … 284 224 DO jk = 1, nlay_s 285 225 DO ji = kideb, kiut 286 zdeltah(ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk) 287 zqfont_su(ji) = MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * & 288 q_s_b(ji,jk) 289 zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) 290 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) !resulting melt of snow 226 zdeltah (ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk) 227 zqfont_su(ji) = MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * q_s_b(ji,jk) 228 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) 229 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) ! resulting melt of snow 291 230 END DO 292 231 END DO … … 302 241 ht_s_b(ji) = MAX( zzero , zhsnew ) 303 242 ! Volume and mass variations of snow 304 dvsbq_1d(ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) & 305 - zdh_s_mel(ji) ) 306 dvsbq_1d(ji) = MIN( zzero, dvsbq_1d(ji) ) 243 dvsbq_1d (ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_mel(ji) ) 244 dvsbq_1d (ji) = MIN( zzero, dvsbq_1d(ji) ) 307 245 rdmsnif_1d(ji) = rdmsnif_1d(ji) + rhosn * dvsbq_1d(ji) 308 246 END DO ! ji … … 312 250 !-------------------------- 313 251 DO ji = kideb, kiut 314 dh_i_surf(ji) = 0.0 315 ! For heat conservation test 316 z_f_surf(ji) = zqfont_su(ji) / rdt_ice ! heat conservation test 317 zdq_i(ji) = 0.0 252 dh_i_surf(ji) = 0.e0 253 z_f_surf (ji) = zqfont_su(ji) / rdt_ice ! heat conservation test 254 zdq_i (ji) = 0.e0 318 255 END DO ! ji 319 256 320 257 DO jk = 1, nlay_i 321 258 DO ji = kideb, kiut 322 ! melt of layer jk 323 zdeltah(ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) 324 ! recompute heat available 325 zqfont_su(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * & 326 q_i_b(ji,jk) 327 ! melt of layer jk cannot be higher than its thickness 328 zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) 329 ! update surface melt 330 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) 331 ! for energy conservation 332 zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * & 333 q_i_b(ji,jk) / rdt_ice 259 ! ! melt of layer jk 260 zdeltah (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) 261 ! ! recompute heat available 262 zqfont_su(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 263 ! ! melt of layer jk cannot be higher than its thickness 264 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) 265 ! ! update surface melt 266 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) 267 ! ! for energy conservation 268 zdq_i (ji) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 269 ! 334 270 ! contribution to ice-ocean salt flux 335 zji = MOD( npb(ji) - 1, jpi ) + 1 336 zjj = ( npb(ji) - 1 ) / jpi + 1 337 zfsalt_melt(ji) = zfsalt_melt(ji) + & 338 ( sss_m(zji,zjj) - sm_i_b(ji) ) * & 339 a_i_b(ji) * & 340 MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 341 END DO ! ji 342 END DO ! jk 343 344 !------------------- 345 ! Conservation test 346 !------------------- 347 IF ( con_i ) THEN 348 numce_dh = 0 349 meance_dh = 0.0 271 zji = MOD( npb(ji) - 1, jpi ) + 1 272 zjj = ( npb(ji) - 1 ) / jpi + 1 273 zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji) ) * a_i_b(ji) & 274 & * MIN( zdeltah(ji,jk) , 0.e0 ) * rhoic / rdt_ice 275 END DO 276 END DO 277 278 ! !------------------- 279 IF( con_i ) THEN ! Conservation test 280 ! !------------------- 281 numce_dh = 0 282 meance_dh = 0.e0 350 283 DO ji = kideb, kiut 351 352 284 IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN 353 285 numce_dh = numce_dh + 1 354 286 meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji) 355 287 ENDIF 356 357 IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN! 288 IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN! 358 289 WRITE(numout,*) ' ALERTE heat loss for surface melt ' 359 290 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl … … 368 299 WRITE(numout,*) ' sss_m : ', sss_m(zji,zjj) 369 300 ENDIF 370 371 END DO ! ji 372 373 IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh 301 END DO 302 ! 303 IF( numce_dh > 0 ) meance_dh = meance_dh / numce_dh 374 304 WRITE(numout,*) ' Error report - Category : ', jl 375 305 WRITE(numout,*) ' ~~~~~~~~~~~~ ' 376 306 WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh 377 307 WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh 378 379 ENDIF ! con_i308 ! 309 ENDIF 380 310 381 311 !---------------------- … … 386 316 ! if qla is positive (upwards), heat goes to the atmosphere, therefore 387 317 ! snow sublimates, if qla is negative (downwards), snow condensates 388 zdh_s_sub(ji) = - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice389 dh_s_tot (ji) = dh_s_tot(ji) + zdh_s_sub(ji)390 zdhcf = ht_s_b(ji) + zdh_s_sub(ji)391 ht_s_b (ji)= MAX( zzero , zdhcf )318 zdh_s_sub(ji) = - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice 319 dh_s_tot (ji) = dh_s_tot(ji) + zdh_s_sub(ji) 320 zdhcf = ht_s_b(ji) + zdh_s_sub(ji) 321 ht_s_b (ji) = MAX( zzero , zdhcf ) 392 322 ! we recompute dh_s_tot 393 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji)394 qt_s_in (ji,jl) = qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1)395 END DO !ji396 397 zqt_dummy(:) = 0. 0323 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 324 qt_s_in (ji,jl) = qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1) 325 END DO 326 327 zqt_dummy(:) = 0.e0 398 328 DO jk = 1, nlay_s 399 329 DO ji = kideb,kiut 400 q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 401 ! heat conservation 402 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s 403 END DO 404 END DO 405 406 DO jk = 1, nlay_s !n 407 DO ji = kideb, kiut !n 408 ! In case of disparition of the snow, we have to update the snow 409 ! temperatures 330 q_s_b (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 331 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s ! heat conservation 332 END DO 333 END DO 334 335 DO jk = 1, nlay_s 336 DO ji = kideb, kiut 337 ! In case of disparition of the snow, we have to update the snow temperatures 410 338 zhisn = MAX( zzero , SIGN( zone, - ht_s_b(ji) ) ) 411 339 t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt … … 428 356 ! 4.1 Basal growth - (a) salinity not varying in time 429 357 !----------------------------------------------------- 430 IF ( ( num_sal .NE. 2 ) .AND. ( num_sal .NE. 4 )) THEN358 IF( num_sal /= 2 .AND. num_sal /= 4 ) THEN 431 359 DO ji = kideb, kiut 432 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT.0.0 ) THEN360 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0.0 ) THEN 433 361 s_i_new(ji) = sm_i_b(ji) 434 362 ! Melting point in K … … 437 365 ztform = t_i_b(ji,nlay_i) ! t_bo_b crashes in the 438 366 ! Baltic 439 q_i_b(ji,nlay_i+1) = rhoic * & 440 ( cpic * ( ztmelts - ztform ) & 441 + lfus * ( 1.0 - ( ztmelts - rtt ) / & 442 ( ztform - rtt ) ) & 443 - rcp * ( ztmelts-rtt ) ) 367 q_i_b(ji,nlay_i+1) = rhoic * ( cpic * ( ztmelts - ztform ) & 368 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( ztform - rtt ) ) & 369 & - rcp * ( ztmelts - rtt ) ) 444 370 ! Basal growth rate = - F*dt / q 445 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & 446 qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 447 ENDIF ! heat budget 448 END DO ! ji 449 ENDIF ! num_sal 371 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 372 ENDIF 373 END DO 374 ENDIF 450 375 451 376 !------------------------------------------------- 452 377 ! 4.1 Basal growth - (b) salinity varying in time 453 378 !------------------------------------------------- 454 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 )) THEN379 IF( num_sal == 2 .OR. num_sal == 4 ) THEN 455 380 ! the growth rate (dh_i_bott) is function of the new ice 456 381 ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice … … 461 386 ! Initial value (tested 1D, can be anything between 1 and 20) 462 387 num_iter_max = 4 463 s_i_new(:) = 4.0388 s_i_new(:) = 4.0 464 389 465 390 ! Iterative procedure 466 391 DO iter = 1, num_iter_max 467 392 DO ji = kideb, kiut 468 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN393 IF( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0 ) THEN 469 394 zji = MOD( npb(ji) - 1, jpi ) + 1 470 395 zjj = ( npb(ji) - 1 ) / jpi + 1 … … 472 397 ztmelts = - tmut * s_i_new(ji) + rtt 473 398 ! New ice heat content (Bitz and Lipscomb, 1999) 474 q_i_b(ji,nlay_i+1) = rhoic * & 475 ( cpic * ( ztmelts - t_bo_b(ji) ) & 476 + lfus * ( 1.0 - ( ztmelts - rtt ) / & 477 ( t_bo_b(ji) - rtt ) ) & 478 - rcp * ( ztmelts-rtt ) ) 399 q_i_b(ji,nlay_i+1) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 400 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) & 401 & - rcp * ( ztmelts-rtt ) ) 479 402 ! Bottom growth rate = - F*dt / q 480 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) & 481 + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 403 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 482 404 ! New ice salinity ( Cox and Weeks, JGR, 1988 ) 483 405 ! zswi2 (1) if dh_i_bott/rdt .GT. 3.6e-7 484 406 ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 485 407 ! zswi1 (1) if dh_i_bott/rdt .LT. 2.0e-8 486 zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , zeps) )408 zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , epsi13 ) ) 487 409 zswi2 = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 488 410 zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 489 411 zswi1 = 1. - zswi2 * zswi12 490 zfracs = zswi1 * 0.12 + & 491 zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + & 492 zswi2 * 0.26 / & 493 ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 494 zds = zfracs*sss_m(zji,zjj) - s_i_new(ji) 412 zfracs = zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) & 413 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 414 zds = zfracs * sss_m(zji,zjj) - s_i_new(ji) 495 415 s_i_new(ji) = zfracs * sss_m(zji,zjj) 496 416 ENDIF ! fc_bo_i … … 500 420 ! Final values 501 421 DO ji = kideb, kiut 502 IF 422 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN 503 423 ! New ice salinity must not exceed 15 psu 504 424 s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) … … 506 426 ztmelts = - tmut * s_i_new(ji) + rtt 507 427 ! New ice heat content (Bitz and Lipscomb, 1999) 508 q_i_b(ji,nlay_i+1) = rhoic * & 509 ( cpic * ( ztmelts - t_bo_b(ji) ) & 510 + lfus * ( 1.0 - ( ztmelts - rtt ) / & 511 ( t_bo_b(ji) - rtt ) ) & 512 - rcp * ( ztmelts-rtt ) ) 428 q_i_b(ji,nlay_i+1) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 429 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) & 430 & - rcp * ( ztmelts - rtt ) ) 513 431 ! Basal growth rate = - F*dt / q 514 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & 515 qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 432 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 516 433 ! Salinity update 517 434 ! entrapment during bottom growth 518 dsm_i_se_1d(ji) = ( s_i_new(ji)*dh_i_bott(ji) + & 519 sm_i_b(ji)*ht_i_b(ji) ) / & 520 MAX( ht_i_b(ji) + dh_i_bott(ji) ,zeps ) & 521 - sm_i_b(ji) 435 dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) ) & 436 & / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji) 522 437 ENDIF ! heat budget 523 END DO ! ji524 ENDIF ! num_sal438 END DO 439 ENDIF 525 440 526 441 !---------------- … … 528 443 !---------------- 529 444 meance_dh = 0.0 530 numce_dh = 0445 numce_dh = 0 531 446 innermelt(:) = 0 532 447 533 448 DO ji = kideb, kiut 534 449 ! heat convergence at the surface > 0 535 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN450 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0 ) THEN 536 451 537 452 s_i_new(ji) = s_i_b(ji,nlay_i) … … 539 454 540 455 zfbase(ji) = zqfont_bo(ji) / rdt_ice ! heat conservation test 541 zdq_i(ji) = 0. 0542 543 dh_i_bott(ji) = 0. 0456 zdq_i(ji) = 0.e0 457 458 dh_i_bott(ji) = 0.e0 544 459 ENDIF 545 460 END DO … … 555 470 ELSE ! normal ablation 556 471 zdeltah(ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) 557 zqfont_bo(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * & 558 q_i_b(ji,jk) 472 zqfont_bo(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 559 473 zdeltah(ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) ) 560 474 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) 561 zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * & 562 q_i_b(ji,jk) / rdt_ice 475 zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 563 476 ! contribution to salt flux 564 477 zji = MOD( npb(ji) - 1, jpi ) + 1 565 478 zjj = ( npb(ji) - 1 ) / jpi + 1 566 zfsalt_melt(ji) = zfsalt_melt(ji) + & 567 ( sss_m(zji,zjj) - sm_i_b(ji) ) * & 568 a_i_b(ji) * & 569 MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 479 zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji) ) * a_i_b(ji) & 480 & * MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 570 481 ENDIF 571 482 ENDIF … … 573 484 END DO ! jk 574 485 575 !------------------- 576 ! Conservation test 577 !------------------- 578 IF ( con_i ) THEN 486 ! !------------------- 487 IF( con_i ) THEN ! Conservation test 488 ! !------------------- 579 489 DO ji = kideb, kiut 580 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN581 IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN582 numce_dh = numce_dh + 1490 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0 ) THEN 491 IF( ( zfbase(ji) + zdq_i(ji) ) >= 1.e-3 ) THEN 492 numce_dh = numce_dh + 1 583 493 meance_dh = meance_dh + zfbase(ji) + zdq_i(ji) 584 494 ENDIF … … 598 508 WRITE(numout,*) ' innermelt : ', innermelt(ji) 599 509 ENDIF 600 ENDIF ! heat convergence at the surface 601 END DO ! ji 602 603 IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh 510 ENDIF 511 END DO 512 IF( numce_dh > 0 ) meance_dh = meance_dh / numce_dh 604 513 WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh 605 514 WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh 606 515 WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d) 607 608 ENDIF ! con_i516 ! 517 ENDIF 609 518 610 519 ! … … 618 527 619 528 DO ji = kideb, kiut 620 ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)621 zdhbf = dh_i_bott(ji)622 IF (jpl.EQ.1) zdhbf = MAX( hmelt , dh_i_bott(ji) )623 ! excessive energy is sent to lateral ablation624 fsup(ji) = rhoic*lfus * at_i_b(ji) / MAX( ( 1.0 - at_i_b(ji) ),epsi13) &625 * ( zdhbf - dh_i_bott(ji) ) / rdt_ice626 529 ! ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) 530 IF( jpl == 1 ) THEN ; zdhbf = MAX( hmelt , dh_i_bott(ji) ) 531 ELSE ; zdhbf = dh_i_bott(ji) 532 ENDIF 533 ! ! excessive energy is sent to lateral ablation 534 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) & 535 & * ( zdhbf - dh_i_bott(ji) ) / rdt_ice 627 536 dh_i_bott(ji) = zdhbf 628 !since ice volume is only used for outputs, we keep it global for all categories 629 dvbbq_1d(ji) = a_i_b(ji)*dh_i_bott(ji) 630 !new ice thickness 631 zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 632 633 ! diagnostic ( bottom ice growth ) 634 zji = MOD( npb(ji) - 1, jpi ) + 1 635 zjj = ( npb(ji) - 1 ) / jpi + 1 636 diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) & 637 / rdt_ice 638 diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) & 639 / rdt_ice 640 diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) & 641 / rdt_ice 537 ! !since ice volume is only used for outputs, we keep it global for all categories 538 dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 539 ! !new ice thickness 540 zhgnew (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 541 ! ! diagnostic ( bottom ice growth ) 542 zji = MOD( npb(ji) - 1, jpi ) + 1 543 zjj = ( npb(ji) - 1 ) / jpi + 1 544 diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice 545 diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) / rdt_ice 546 diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice 642 547 END DO 643 548 … … 653 558 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 654 559 ! 0 if no more ice 655 zhgnew (ji) = zihgnew * zhgnew(ji)! ice thickness is put to 0560 zhgnew (ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 656 561 ! remaining heat 657 562 zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) 658 563 659 564 ! If snow remains, energy is used to melt snow 660 zhni = ht_s_b(ji)! snow depth at previous time step661 zihg 565 zhni = ht_s_b(ji) ! snow depth at previous time step 566 zihg = MAX( zzero , SIGN ( zone , - ht_s_b(ji) ) ) ! 0 if snow 662 567 663 568 ! energy of melting of remaining snow 664 zqt_s(ji) = ( 1. - zihg) * zqt_s(ji) / MAX( zhni, zeps ) 665 zdhnm = - ( 1. - zihg ) * ( 1. - zihgnew ) * ( zfdt_final(ji) / & 666 MAX( zqt_s(ji) , zeps ) ) 569 zqt_s(ji) = ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi13 ) 570 zdhnm = - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 ) 667 571 zhnfi = zhni + zdhnm 668 zfdt_final(ji) = MAX 572 zfdt_final(ji) = MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 ) 669 573 ht_s_b(ji) = MAX( zzero , zhnfi ) 670 574 zqt_s(ji) = zqt_s(ji) * ht_s_b(ji) … … 672 576 ! Mass variations of ice and snow 673 577 !--------------------------------- 674 ! 578 ! ! mass variation of the jl category 675 579 zzfmass_s = - a_i_b(ji) * ( zhni - ht_s_b(ji) ) * rhosn ! snow 676 580 zzfmass_i = a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic ! ice … … 684 588 ! Remaining heat to the ocean 685 589 !--------------------------------- 686 ! focea is in W.m-2 * dt 687 focea(ji) = - zfdt_final(ji) / rdt_ice 590 focea(ji) = - zfdt_final(ji) / rdt_ice ! focea is in W.m-2 * dt 688 591 689 592 END DO … … 695 598 !--------------------------- 696 599 DO ji = kideb, kiut 697 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice600 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 698 601 699 602 ! Salt flux 700 zji = MOD( npb(ji) - 1, jpi ) + 1 701 zjj = ( npb(ji) - 1 ) / jpi + 1 702 IF ( num_sal .NE. 4 ) & 703 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + & 704 (1.0 - zihgnew) * zfmass_i(ji) * & 705 ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice 603 zji = MOD( npb(ji) - 1, jpi ) + 1 604 zjj = ( npb(ji) - 1 ) / jpi + 1 706 605 ! new lines 707 IF ( num_sal .EQ. 4 ) & 708 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + & 709 (1.0 - zihgnew) * zfmass_i(ji) * & 710 ( sss_m(zji,zjj) - bulk_sal ) / rdt_ice 606 IF( num_sal == 4 ) THEN 607 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) & 608 & + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - bulk_sal ) / rdt_ice 609 ELSE 610 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) & 611 & + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice 612 ENDIF 711 613 ! Heat flux 712 614 ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1 713 615 ! excessive total ablation energy (focea) sent to the ocean 714 qfvbq_1d(ji) = qfvbq_1d(ji) + & 715 fsup(ji) + ( 1.0 - zihgnew ) * & 716 focea(ji) * a_i_b(ji) * rdt_ice 616 qfvbq_1d(ji) = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice 717 617 718 618 zihic = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) ) 719 619 ! equals 0 if ht_i = 0, 1 if ht_i gt 0 720 620 fscbq_1d(ji) = a_i_b(ji) * fstbif_1d(ji) 721 qldif_1d(ji) = qldif_1d(ji) & 722 + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) & 723 * rdt_ice & 724 + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice 621 qldif_1d(ji) = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice & 622 & + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice 725 623 END DO ! ji 726 624 … … 729 627 !------------------------------------------- 730 628 DO ji = kideb, kiut 731 zihgnew 732 t_su_b(ji) 629 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 630 t_su_b(ji) = zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt 733 631 END DO ! ji 734 632 735 633 DO jk = 1, nlay_i 736 634 DO ji = kideb, kiut 737 zihgnew 738 t_i_b(ji,jk) 739 q_i_b(ji,jk) 635 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 636 t_i_b(ji,jk) = zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt 637 q_i_b(ji,jk) = zihgnew * q_i_b(ji,jk) 740 638 END DO 741 639 END DO ! ji … … 748 646 ! 6) Snow-Ice formation | 749 647 !------------------------------------------------------------------------------| 750 ! 751 ! When snow load excesses Archimede's limit, snow-ice interface goes down 752 ! under sea-level, flooding of seawater transforms snow into ice 753 ! dh_snowice is positive for the ice 754 DO ji = kideb, kiut 755 756 dh_snowice(ji) = MAX(zzero,(rhosn*ht_s_b(ji)+(rhoic-rau0) & 757 * ht_i_b(ji))/(rhosn+rau0-rhoic)) 758 zhgnew(ji) = MAX(zhgnew(ji),zhgnew(ji)+dh_snowice(ji)) 759 zhnnew = MIN(ht_s_b(ji),ht_s_b(ji)-dh_snowice(ji)) 648 ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level, 649 ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 650 DO ji = kideb, kiut 651 ! 652 dh_snowice(ji) = MAX( zzero , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic ) ) 653 zhgnew(ji) = MAX( zhgnew(ji) , zhgnew(ji) + dh_snowice(ji) ) 654 zhnnew = MIN( ht_s_b(ji) , ht_s_b(ji) - dh_snowice(ji) ) 760 655 761 656 ! Changes in ice volume and ice mass. 762 dvnbq_1d(ji) = a_i_b(ji) * (zhgnew(ji)-ht_i_b(ji)) 763 dmgwi_1d(ji) = dmgwi_1d(ji) + a_i_b(ji) & 764 *(ht_s_b(ji)-zhnnew)*rhosn 765 766 rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) & 767 * ( zhgnew(ji) - ht_i_b(ji) )*rhoic 768 rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) & 769 * ( zhnnew - ht_s_b(ji) )*rhosn 657 dvnbq_1d (ji) = a_i_b(ji) * ( zhgnew(ji)-ht_i_b(ji) ) 658 dmgwi_1d (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 659 660 rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic 661 rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn 770 662 771 663 ! Equivalent salt flux (1) Snow-ice formation component 772 664 ! ----------------------------------------------------- 773 zji = MOD( npb(ji) - 1, jpi ) + 1 774 zjj = ( npb(ji) - 1 ) / jpi + 1 775 776 zsm_snowice = ( rhoic - rhosn ) / rhoic * & 777 sss_m(zji,zjj) 778 779 IF ( num_sal .NE. 2 ) zsm_snowice = sm_i_b(ji) 780 781 IF ( num_sal .NE. 4 ) & 782 fseqv_1d(ji) = fseqv_1d(ji) + & 783 ( sss_m(zji,zjj) - zsm_snowice ) * & 784 a_i_b(ji) * & 785 ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 786 ! new lines 787 IF ( num_sal .EQ. 4 ) & 788 fseqv_1d(ji) = fseqv_1d(ji) + & 789 ( sss_m(zji,zjj) - bulk_sal ) * & 790 a_i_b(ji) * & 791 ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 792 665 zji = MOD( npb(ji) - 1, jpi ) + 1 666 zjj = ( npb(ji) - 1 ) / jpi + 1 667 668 IF( num_sal /= 2 ) THEN ; zsm_snowice = sm_i_b(ji) 669 ELSE ; zsm_snowice = ( rhoic - rhosn ) / rhoic * sss_m(zji,zjj) 670 ENDIF 671 IF( num_sal == 4 ) THEN 672 fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal ) * a_i_b(ji) & 673 & * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 674 ELSE 675 fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - zsm_snowice ) * a_i_b(ji) & 676 & * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 677 ENDIF 793 678 ! entrapment during snow ice formation 794 i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 795 isnowic = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * & 796 i_ice_switch 797 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 798 dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) & 799 + sm_i_b(ji) * ht_i_b(ji) & 800 / MAX( ht_i_b(ji) + dh_snowice(ji), zeps) & 801 - sm_i_b(ji) ) * isnowic 679 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 680 isnowic = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji) ) ) * i_ice_switch 681 IF( num_sal == 2 .OR. num_sal == 4 ) & 682 dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) & 683 & + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13) & 684 & - sm_i_b(ji) ) * isnowic 802 685 803 686 ! Actualize new snow and ice thickness. … … 806 689 807 690 ! Total ablation ! new lines added to debug 808 IF( ht_i_b(ji) .LE.0.0 )a_i_b(ji) = 0.0691 IF( ht_i_b(ji) <= 0.e0 ) a_i_b(ji) = 0.0 809 692 810 693 ! diagnostic ( snow ice growth ) 811 zji = MOD( npb(ji) - 1, jpi ) + 1 812 zjj = ( npb(ji) - 1 ) / jpi + 1 813 diag_sni_gr(zji,zjj) = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / & 814 rdt_ice 815 694 zji = MOD( npb(ji) - 1, jpi ) + 1 695 zjj = ( npb(ji) - 1 ) / jpi + 1 696 diag_sni_gr(zji,zjj) = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / rdt_ice 697 ! 816 698 END DO !ji 817 699 818 700 END SUBROUTINE lim_thd_dh 701 819 702 #else 820 !!====================================================================== 821 !! *** MODULE limthd_dh *** 822 !! no sea ice model 823 !!====================================================================== 703 !!---------------------------------------------------------------------- 704 !! Default option NO LIM3 sea-ice model 705 !!---------------------------------------------------------------------- 824 706 CONTAINS 825 707 SUBROUTINE lim_thd_dh ! Empty routine 826 708 END SUBROUTINE lim_thd_dh 827 709 #endif 710 711 !!====================================================================== 828 712 END MODULE limthd_dh
Note: See TracChangeset
for help on using the changeset viewer.
