Changeset 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
- Timestamp:
- 2015-12-01T16:35:30+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4624 r5965 8 8 !! 3.0 ! 2005-11 (M. Vancoppenolle) LIM-3 : Multi-layer thermodynamics + salinity variations 9 9 !! - ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec, lim_thd_con_dh and lim_thd_con_dif 10 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_snw10 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw 11 11 !! 3.3 ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas) 12 12 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation … … 22 22 USE phycst ! physical constants 23 23 USE dom_oce ! ocean space and time domain variables 24 USE oce , ONLY : iatte, oatte25 24 USE ice ! LIM: sea-ice variables 26 USE par_ice ! LIM: sea-ice parameters27 25 USE sbc_oce ! Surface boundary condition: ocean fields 28 26 USE sbc_ice ! Surface boundary condition: ice fields 29 27 USE thd_ice ! LIM thermodynamic sea-ice variables 30 28 USE dom_ice ! LIM sea-ice domain 31 USE domvvl ! domain: variable volume level32 29 USE limthd_dif ! LIM: thermodynamics, vertical diffusion 33 30 USE limthd_dh ! LIM: thermodynamics, ice and snow thickness variation 34 31 USE limthd_sal ! LIM: thermodynamics, ice salinity 35 32 USE limthd_ent ! LIM: thermodynamics, ice enthalpy redistribution 33 USE limthd_lac ! LIM-3 lateral accretion 34 USE limitd_th ! remapping thickness distribution 36 35 USE limtab ! LIM: 1D <==> 2D transformation 37 36 USE limvar ! LIM: sea-ice variables … … 43 42 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 44 43 USE timing ! Timing 44 USE limcons ! conservation tests 45 USE limctl 45 46 46 47 IMPLICIT NONE 47 48 PRIVATE 48 49 49 PUBLIC lim_thd ! called by limstp module 50 PUBLIC lim_thd_init ! called by iceini module 51 52 REAL(wp) :: epsi10 = 1.e-10_wp ! 53 REAL(wp) :: zzero = 0._wp ! 54 REAL(wp) :: zone = 1._wp ! 50 PUBLIC lim_thd ! called by limstp module 51 PUBLIC lim_thd_init ! called by sbc_lim_init 55 52 56 53 !! * Substitutions … … 68 65 !! *** ROUTINE lim_thd *** 69 66 !! 70 !! ** Purpose : This routine manages the ice thermodynamic.67 !! ** Purpose : This routine manages ice thermodynamics 71 68 !! 72 69 !! ** Action : - Initialisation of some variables … … 74 71 !! at the ice base, snow acc.,heat budget of the leads) 75 72 !! - selection of the icy points and put them in an array 76 !! - call lim_vert_ther for vert ice thermodynamic 77 !! - back to the geographic grid 78 !! - selection of points for lateral accretion 79 !! - call lim_lat_acc for the ice accretion 73 !! - call lim_thd_dif for vertical heat diffusion 74 !! - call lim_thd_dh for vertical ice growth and melt 75 !! - call lim_thd_ent for enthalpy remapping 76 !! - call lim_thd_sal for ice desalination 77 !! - call lim_thd_temp to retrieve temperature from ice enthalpy 80 78 !! - back to the geographic grid 81 79 !! 82 !! ** References : H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-9080 !! ** References : 83 81 !!--------------------------------------------------------------------- 84 INTEGER, INTENT(in) :: 82 INTEGER, INTENT(in) :: kt ! number of iteration 85 83 !! 86 INTEGER :: ji, jj, jk, jl ! dummy loop indices 87 INTEGER :: nbpb ! nb of icy pts for thermo. cal. 88 REAL(wp) :: zfric_umin = 5e-03_wp ! lower bound for the friction velocity 89 REAL(wp) :: zfric_umax = 2e-02_wp ! upper bound for the friction velocity 90 REAL(wp) :: zinda, zindb, zthsnice, zfric_u ! local scalar 91 REAL(wp) :: zfntlat, zpareff, zareamin, zcoef ! - - 92 REAL(wp), POINTER, DIMENSION(:,:) :: zqlbsbq ! link with lead energy budget qldif 93 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 94 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 84 INTEGER :: ji, jj, jk, jl ! dummy loop indices 85 INTEGER :: nbpb ! nb of icy pts for vertical thermo calculations 86 INTEGER :: ii, ij ! temporary dummy loop index 87 REAL(wp) :: zfric_u, zqld, zqfr 88 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 89 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 90 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 91 ! 95 92 !!------------------------------------------------------------------- 93 96 94 IF( nn_timing == 1 ) CALL timing_start('limthd') 97 95 98 CALL wrk_alloc( jpi, jpj, zqlbsbq ) 99 100 ! ------------------------------- 101 !- check conservation (C Rousset) 102 IF (ln_limdiahsb) THEN 103 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 104 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 105 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 106 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 107 ENDIF 108 !- check conservation (C Rousset) 109 ! ------------------------------- 110 111 !------------------------------------------------------------------------------! 112 ! 1) Initialization of diagnostic variables ! 113 !------------------------------------------------------------------------------! 96 ! conservation test 97 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 98 99 CALL lim_var_glo2eqv 100 !------------------------------------------------------------------------! 101 ! 1) Initialization of some variables ! 102 !------------------------------------------------------------------------! 103 ftr_ice(:,:,:) = 0._wp ! part of solar radiation transmitted through the ice 114 104 115 105 !-------------------- 116 106 ! 1.2) Heat content 117 107 !-------------------- 118 ! Change the units of heat content; from global units to J.m3108 ! Change the units of heat content; from J/m2 to J/m3 119 109 DO jl = 1, jpl 120 110 DO jk = 1, nlay_i 121 111 DO jj = 1, jpj 122 112 DO ji = 1, jpi 113 !0 if no ice and 1 if yes 114 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 ) ) 123 115 !Energy of melting q(S,T) [J.m-3] 124 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 125 !0 if no ice and 1 if yes 126 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) 127 !convert units ! very important that this line is here 128 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 116 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i ) 129 117 END DO 130 118 END DO … … 133 121 DO jj = 1, jpj 134 122 DO ji = 1, jpi 123 !0 if no ice and 1 if yes 124 rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) ) 135 125 !Energy of melting q(S,T) [J.m-3] 136 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s ) 137 !0 if no ice and 1 if yes 138 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 ) ) 139 !convert units ! very important that this line is here 140 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb 126 e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s ) 141 127 END DO 142 128 END DO … … 144 130 END DO 145 131 146 !-----------------------------------147 ! 1.4) Compute global heat content148 !-----------------------------------149 qt_i_in (:,:) = 0.e0150 qt_s_in (:,:) = 0.e0151 qt_i_fin (:,:) = 0.e0152 qt_s_fin (:,:) = 0.e0153 sum_fluxq(:,:) = 0.e0154 fatm (:,:) = 0.e0155 156 132 ! 2) Partial computation of forcing for the thermodynamic sea ice model. ! 157 133 !-----------------------------------------------------------------------------! 158 159 !CDIR NOVERRCHK160 134 DO jj = 1, jpj 161 !CDIR NOVERRCHK162 135 DO ji = 1, jpi 163 zinda = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) + epsi10 ) ) )136 rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 164 137 ! 165 138 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 168 141 ! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean 169 142 ! ! temperature and turbulent mixing (McPhee, 1992) 170 ! friction velocity171 zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )172 173 ! here the drag will depend on ice thickness and type (0.006)174 fdtcn(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) )175 ! also category dependent176 ! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead177 qdtcn(ji,jj) = zinda * fdtcn(ji,jj) * ( 1.0 - at_i(ji,jj) ) * rdt_ice178 !179 ! !-- Lead heat budget, qldif (part 1, next one is in limthd_dh)180 ! ! caution: exponent betas used as more snow can fallinto leads181 qldif(ji,jj) = tms(ji,jj) * rdt_ice * ( &182 & pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif183 & + qns(ji,jj) & ! non solar heat184 & + fdtcn(ji,jj) & ! turbulent ice-ocean heat185 & + fsbbq(ji,jj) * ( 1.0 - zinda ) ) & ! residual heat from previous step186 & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus ) ! latent heat of sprecip melting187 143 ! 188 ! Positive heat budget is used for bottom ablation 189 zfntlat = 1.0 - MAX( zzero , SIGN( zone , - qldif(ji,jj) ) ) 190 != 1 if positive heat budget 191 zpareff = 1.0 - zinda * zfntlat 192 != 0 if ice and positive heat budget and 1 if one of those two is false 193 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / ( rdt_ice * MAX( at_i(ji,jj), epsi10 ) ) 144 ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- ! 145 zqld = tmask(ji,jj,1) * rdt_ice * & 146 & ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 147 148 ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 149 zqfr = tmask(ji,jj,1) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 150 151 ! --- Energy from the turbulent oceanic heat flux (W/m2) --- ! 152 zfric_u = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 153 fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2 154 fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 155 ! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach 156 ! the freezing point, so that we do not have SST < T_freeze 157 ! This implies: - ( fhtur(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0 158 159 !-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice 160 qlead(ji,jj) = MIN( 0._wp , zqld - ( fhtur(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr ) 161 162 ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting 163 IF( zqld > 0._wp ) THEN 164 fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90 165 qlead(ji,jj) = 0._wp 166 ELSE 167 fhld (ji,jj) = 0._wp 168 ENDIF 194 169 ! 195 ! Heat budget of the lead, energy transferred from ice to ocean 196 qldif (ji,jj) = zpareff * qldif(ji,jj) 197 qdtcn (ji,jj) = zpareff * qdtcn(ji,jj) 198 ! 199 ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 200 qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 201 ! 202 ! oceanic heat flux (limthd_dh) 203 fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) + fdtcn(ji,jj) ) 204 ! 170 ! ----------------------------------------- 171 ! Net heat flux on top of ice-ocean [W.m-2] 172 ! ----------------------------------------- 173 hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 174 175 ! ----------------------------------------------------------------------------- 176 ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2] 177 ! ----------------------------------------------------------------------------- 178 ! First step here : non solar + precip - qlead - qturb 179 ! Second step in limthd_dh : heat remaining if total melt (zq_rema) 180 ! Third step in limsbc : heat from ice-ocean mass exchange (zf_mass) + solar 181 hfx_out(ji,jj) = pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) & ! Non solar heat flux received by the ocean 182 & - qlead(ji,jj) * r1_rdtice & ! heat flux taken from the ocean where there is open water ice formation 183 & - at_i(ji,jj) * fhtur(ji,jj) & ! heat flux taken by turbulence 184 & - at_i(ji,jj) * fhld(ji,jj) ! heat flux taken during bottom growth/melt 185 ! (fhld should be 0 while bott growth) 205 186 END DO 206 187 END DO … … 217 198 ENDIF 218 199 219 zareamin = epsi10220 200 nbpb = 0 221 201 DO jj = 1, jpj 222 202 DO ji = 1, jpi 223 IF ( a_i(ji,jj,jl) .gt. zareamin) THEN203 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 224 204 nbpb = nbpb + 1 225 205 npb(nbpb) = (jj - 1) * jpi + ji … … 230 210 ! debug point to follow 231 211 jiindex_1d = 0 232 IF( ln_ nicep) THEN233 DO ji = mi0( jiindx), mi1(jiindx)234 DO jj = mj0(j jindx), mj1(jjindx)212 IF( ln_icectl ) THEN 213 DO ji = mi0(iiceprt), mi1(iiceprt) 214 DO jj = mj0(jiceprt), mj1(jiceprt) 235 215 jiindex_1d = (jj - 1) * jpi + ji 216 WRITE(numout,*) ' lim_thd : Category no : ', jl 236 217 END DO 237 218 END DO … … 246 227 IF( nbpb > 0 ) THEN ! If there is no ice, do nothing. 247 228 248 !------------------------- 249 ! 4.1 Move to 1D arrays 250 !------------------------- 251 252 CALL tab_2d_1d( nbpb, at_i_b (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) ) 253 CALL tab_2d_1d( nbpb, a_i_b (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 254 CALL tab_2d_1d( nbpb, ht_i_b (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 255 CALL tab_2d_1d( nbpb, ht_s_b (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 256 257 CALL tab_2d_1d( nbpb, t_su_b (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 258 CALL tab_2d_1d( nbpb, sm_i_b (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 259 DO jk = 1, nlay_s 260 CALL tab_2d_1d( nbpb, t_s_b(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 261 CALL tab_2d_1d( nbpb, q_s_b(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 262 END DO 263 DO jk = 1, nlay_i 264 CALL tab_2d_1d( nbpb, t_i_b(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 265 CALL tab_2d_1d( nbpb, q_i_b(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 266 CALL tab_2d_1d( nbpb, s_i_b(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 267 END DO 268 269 CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:) , jpi, jpj, npb(1:nbpb) ) 270 CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 271 CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) ) 272 CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb), fr2_i0 , jpi, jpj, npb(1:nbpb) ) 273 CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 274 #if ! defined key_coupled 275 CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 276 CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 277 #endif 278 CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 279 CALL tab_2d_1d( nbpb, t_bo_b (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) ) 280 CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) ) 281 CALL tab_2d_1d( nbpb, fbif_1d (1:nbpb), fbif , jpi, jpj, npb(1:nbpb) ) 282 CALL tab_2d_1d( nbpb, qldif_1d (1:nbpb), qldif , jpi, jpj, npb(1:nbpb) ) 283 CALL tab_2d_1d( nbpb, rdm_ice_1d (1:nbpb), rdm_ice , jpi, jpj, npb(1:nbpb) ) 284 CALL tab_2d_1d( nbpb, rdm_snw_1d (1:nbpb), rdm_snw , jpi, jpj, npb(1:nbpb) ) 285 CALL tab_2d_1d( nbpb, dmgwi_1d (1:nbpb), dmgwi , jpi, jpj, npb(1:nbpb) ) 286 CALL tab_2d_1d( nbpb, qlbbq_1d (1:nbpb), zqlbsbq , jpi, jpj, npb(1:nbpb) ) 287 288 CALL tab_2d_1d( nbpb, sfx_thd_1d (1:nbpb), sfx_thd , jpi, jpj, npb(1:nbpb) ) 289 CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) ) 290 CALL tab_2d_1d( nbpb, fhbri_1d (1:nbpb), fhbri , jpi, jpj, npb(1:nbpb) ) 291 CALL tab_2d_1d( nbpb, fstbif_1d (1:nbpb), fstric , jpi, jpj, npb(1:nbpb) ) 292 CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb), qfvbq , jpi, jpj, npb(1:nbpb) ) 293 294 CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) ) ! clem modif 295 CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte , jpi, jpj, npb(1:nbpb) ) ! clem modif 296 !-------------------------------- 297 ! 4.3) Thermodynamic processes 298 !-------------------------------- 299 300 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting 301 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl ) 302 303 ! !---------------------------------! 304 CALL lim_thd_dif( 1, nbpb, jl ) ! Ice/Snow Temperature profile ! 305 ! !---------------------------------! 306 307 CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting compulsory for limthd_dh 308 309 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 310 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_con_dif( 1 , nbpb , jl ) 311 312 ! !---------------------------------! 313 CALL lim_thd_dh( 1, nbpb, jl ) ! Ice/Snow thickness ! 314 ! !---------------------------------! 315 316 ! !---------------------------------! 317 CALL lim_thd_ent( 1, nbpb, jl ) ! Ice/Snow enthalpy remapping ! 318 ! !---------------------------------! 319 320 ! !---------------------------------! 321 CALL lim_thd_sal( 1, nbpb ) ! Ice salinity computation ! 322 ! !---------------------------------! 323 324 ! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 325 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 326 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_con_dh ( 1 , nbpb , jl ) 327 328 !-------------------------------- 329 ! 4.4) Move 1D to 2D vectors 330 !-------------------------------- 331 332 CALL tab_1d_2d( nbpb, at_i , npb, at_i_b (1:nbpb) , jpi, jpj ) 333 CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_b (1:nbpb) , jpi, jpj ) 334 CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_b (1:nbpb) , jpi, jpj ) 335 CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_b (1:nbpb) , jpi, jpj ) 336 CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_b (1:nbpb) , jpi, jpj ) 337 CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_b (1:nbpb) , jpi, jpj ) 338 DO jk = 1, nlay_s 339 CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b (1:nbpb,jk), jpi, jpj) 340 CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_b (1:nbpb,jk), jpi, jpj) 341 END DO 342 DO jk = 1, nlay_i 343 CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b (1:nbpb,jk), jpi, jpj) 344 CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b (1:nbpb,jk), jpi, jpj) 345 CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b (1:nbpb,jk), jpi, jpj) 346 END DO 347 CALL tab_1d_2d( nbpb, fstric , npb, fstbif_1d (1:nbpb) , jpi, jpj ) 348 CALL tab_1d_2d( nbpb, qldif , npb, qldif_1d (1:nbpb) , jpi, jpj ) 349 CALL tab_1d_2d( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb) , jpi, jpj ) 350 CALL tab_1d_2d( nbpb, rdm_ice , npb, rdm_ice_1d(1:nbpb) , jpi, jpj ) 351 CALL tab_1d_2d( nbpb, rdm_snw , npb, rdm_snw_1d(1:nbpb) , jpi, jpj ) 352 CALL tab_1d_2d( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb) , jpi, jpj ) 353 CALL tab_1d_2d( nbpb, rdvosif , npb, dvsbq_1d (1:nbpb) , jpi, jpj ) 354 CALL tab_1d_2d( nbpb, rdvobif , npb, dvbbq_1d (1:nbpb) , jpi, jpj ) 355 CALL tab_1d_2d( nbpb, fdvolif , npb, dvlbq_1d (1:nbpb) , jpi, jpj ) 356 CALL tab_1d_2d( nbpb, rdvonif , npb, dvnbq_1d (1:nbpb) , jpi, jpj ) 357 CALL tab_1d_2d( nbpb, sfx_thd , npb, sfx_thd_1d(1:nbpb) , jpi, jpj ) 358 ! 359 IF( num_sal == 2 ) THEN 360 CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj ) 361 CALL tab_1d_2d( nbpb, fhbri , npb, fhbri_1d (1:nbpb) , jpi, jpj ) 362 ENDIF 363 ! 364 !+++++ temporary stuff for a dummy version 365 CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb) , jpi, jpj ) 366 CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb) , jpi, jpj ) 367 CALL tab_1d_2d( nbpb, fsup2D , npb, fsup (1:nbpb) , jpi, jpj ) 368 CALL tab_1d_2d( nbpb, focea2D , npb, focea (1:nbpb) , jpi, jpj ) 369 CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new (1:nbpb) , jpi, jpj ) 370 CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0 (1:nbpb) , jpi, jpj ) 371 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qnsr_ice_1d(1:nbpb), jpi, jpj) 372 !+++++ 229 !-------------------------! 230 ! --- Move to 1D arrays --- 231 !-------------------------! 232 CALL lim_thd_1d2d( nbpb, jl, 1 ) 233 234 !--------------------------------------! 235 ! --- Ice/Snow Temperature profile --- ! 236 !--------------------------------------! 237 CALL lim_thd_dif( 1, nbpb ) 238 239 !---------------------------------! 240 ! --- Ice/Snow thickness --- ! 241 !---------------------------------! 242 CALL lim_thd_dh( 1, nbpb ) 243 244 ! --- Ice enthalpy remapping --- ! 245 CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) 246 247 !---------------------------------! 248 ! --- Ice salinity --- ! 249 !---------------------------------! 250 CALL lim_thd_sal( 1, nbpb ) 251 252 !---------------------------------! 253 ! --- temperature update --- ! 254 !---------------------------------! 255 CALL lim_thd_temp( 1, nbpb ) 256 257 !------------------------------------! 258 ! --- lateral melting if monocat --- ! 259 !------------------------------------! 260 IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN 261 CALL lim_thd_lam( 1, nbpb ) 262 END IF 263 264 !-------------------------! 265 ! --- Move to 2D arrays --- 266 !-------------------------! 267 CALL lim_thd_1d2d( nbpb, jl, 2 ) 268 373 269 ! 374 270 IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? 375 271 ENDIF 376 272 ! 377 END DO 273 END DO !jl 378 274 379 275 !------------------------------------------------------------------------------! … … 382 278 383 279 !------------------------ 384 ! 5.1)Ice heat content280 ! Ice heat content 385 281 !------------------------ 386 ! Enthalpies are global variables we have to readjust the units (heat content in 10^9 Joules) 387 zcoef = 1._wp / ( unit_fac * REAL( nlay_i ) ) 282 ! Enthalpies are global variables we have to readjust the units (heat content in J/m2) 388 283 DO jl = 1, jpl 389 284 DO jk = 1, nlay_i 390 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a rea(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) * zcoef285 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i 391 286 END DO 392 287 END DO 393 288 394 289 !------------------------ 395 ! 5.2)Snow heat content290 ! Snow heat content 396 291 !------------------------ 397 ! Enthalpies are global variables we have to readjust the units (heat content in 10^9 Joules) 398 zcoef = 1._wp / ( unit_fac * REAL( nlay_s ) ) 292 ! Enthalpies are global variables we have to readjust the units (heat content in J/m2) 399 293 DO jl = 1, jpl 400 294 DO jk = 1, nlay_s 401 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a rea(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) * zcoef295 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s 402 296 END DO 403 297 END DO 404 298 405 299 !---------------------------------- 406 ! 5.3)Change thickness to volume300 ! Change thickness to volume 407 301 !---------------------------------- 408 CALL lim_var_eqv2glo 302 v_i(:,:,:) = ht_i(:,:,:) * a_i(:,:,:) 303 v_s(:,:,:) = ht_s(:,:,:) * a_i(:,:,:) 304 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 305 306 ! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting in monocat) 307 DO jl = 1, jpl 308 DO jj = 1, jpj 309 DO ji = 1, jpi 310 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi10 ) ) 311 oa_i(ji,jj,jl) = rswitch * oa_i(ji,jj,jl) * a_i(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) 312 END DO 313 END DO 314 END DO 315 316 CALL lim_var_zapsmall 409 317 410 318 !-------------------------------------------- 411 ! 5.4)Diagnostic thermodynamic growth rates319 ! Diagnostic thermodynamic growth rates 412 320 !-------------------------------------------- 413 !clem@useless d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes 414 !clem@mv-to-itd dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 415 416 IF( con_i .AND. jiindex_1d > 0 ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 321 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' ) ! control print 417 322 418 323 IF(ln_ctl) THEN ! Control print … … 420 325 CALL prt_ctl_info(' - Cell values : ') 421 326 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 422 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_thd : cell area :')327 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_thd : cell area :') 423 328 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd : at_i :') 424 329 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd : vt_i :') … … 448 353 ENDIF 449 354 ! 450 ! ------------------------------- 451 !- check conservation (C Rousset) 452 IF (ln_limdiahsb) THEN 453 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 454 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 355 ! 356 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 357 358 !------------------------------------------------------------------------------| 359 ! 6) Transport of ice between thickness categories. | 360 !------------------------------------------------------------------------------| 361 ! Given thermodynamic growth rates, transport ice between thickness categories. 362 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 363 364 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt ) 365 366 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 367 368 !------------------------------------------------------------------------------| 369 ! 7) Add frazil ice growing in leads. 370 !------------------------------------------------------------------------------| 371 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 372 373 CALL lim_thd_lac 374 375 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 376 377 ! Control print 378 IF(ln_ctl) THEN 379 CALL lim_var_glo2eqv 380 381 CALL prt_ctl_info(' ') 382 CALL prt_ctl_info(' - Cell values : ') 383 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 384 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_th : cell area :') 385 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th : at_i :') 386 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th : vt_i :') 387 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th : vt_s :') 388 DO jl = 1, jpl 389 CALL prt_ctl_info(' ') 390 CALL prt_ctl_info(' - Category : ', ivar1=jl) 391 CALL prt_ctl_info(' ~~~~~~~~~~') 392 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_th : a_i : ') 393 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_th : ht_i : ') 394 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_th : ht_s : ') 395 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_th : v_i : ') 396 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_th : v_s : ') 397 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_th : e_s : ') 398 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_th : t_su : ') 399 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_th : t_snow : ') 400 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_th : sm_i : ') 401 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_th : smv_i : ') 402 DO jk = 1, nlay_i 403 CALL prt_ctl_info(' ') 404 CALL prt_ctl_info(' - Layer : ', ivar1=jk) 405 CALL prt_ctl_info(' ~~~~~~~') 406 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : t_i : ') 407 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : e_i : ') 408 END DO 409 END DO 410 ENDIF 411 ! 412 IF( nn_timing == 1 ) CALL timing_stop('limthd') 413 414 END SUBROUTINE lim_thd 415 455 416 456 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 457 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 458 459 zchk_vmin = glob_min(v_i) 460 zchk_amax = glob_max(SUM(a_i,dim=3)) 461 zchk_amin = glob_min(a_i) 462 463 IF(lwp) THEN 464 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limthd) = ',(zchk_v_i * rday) 465 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * rday) 466 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limthd) = ',(zchk_vmin * 1.e-3) 467 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limthd) = ',zchk_amax 468 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limthd) = ',zchk_amin 469 ENDIF 470 ENDIF 471 !- check conservation (C Rousset) 472 ! ------------------------------- 473 ! 474 CALL wrk_dealloc( jpi, jpj, zqlbsbq ) 475 ! 476 IF( nn_timing == 1 ) CALL timing_stop('limthd') 477 END SUBROUTINE lim_thd 478 479 480 SUBROUTINE lim_thd_glohec( eti, ets, etilayer, kideb, kiut, jl ) 417 SUBROUTINE lim_thd_temp( kideb, kiut ) 481 418 !!----------------------------------------------------------------------- 482 !! *** ROUTINE lim_thd_ glohec***419 !! *** ROUTINE lim_thd_temp *** 483 420 !! 484 !! ** Purpose : Compute total heat content for each category 485 !! Works with 1d vectors only 486 !!----------------------------------------------------------------------- 487 INTEGER , INTENT(in ) :: kideb, kiut ! bounds for the spatial loop 488 INTEGER , INTENT(in ) :: jl ! category number 489 REAL(wp), INTENT( out), DIMENSION (jpij,jpl ) :: eti, ets ! vertically-summed heat content for ice & snow 490 REAL(wp), INTENT( out), DIMENSION (jpij,jkmax) :: etilayer ! heat content for ice layers 491 !! 492 INTEGER :: ji,jk ! loop indices 493 !!----------------------------------------------------------------------- 494 eti(:,:) = 0._wp 495 ets(:,:) = 0._wp 496 ! 497 DO jk = 1, nlay_i ! total q over all layers, ice [J.m-2] 498 DO ji = kideb, kiut 499 etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 500 eti (ji,jl) = eti(ji,jl) + etilayer(ji,jk) 501 END DO 502 END DO 503 DO ji = kideb, kiut ! total q over all layers, snow [J.m-2] 504 ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / REAL( nlay_s ) 505 END DO 506 ! 507 WRITE(numout,*) ' lim_thd_glohec ' 508 WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice 509 WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice 510 WRITE(numout,*) ' qt_in : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice 511 ! 512 END SUBROUTINE lim_thd_glohec 513 514 515 SUBROUTINE lim_thd_con_dif( kideb, kiut, jl ) 516 !!----------------------------------------------------------------------- 517 !! *** ROUTINE lim_thd_con_dif *** 518 !! 519 !! ** Purpose : Test energy conservation after heat diffusion 520 !!------------------------------------------------------------------- 521 INTEGER , INTENT(in ) :: kideb, kiut ! bounds for the spatial loop 522 INTEGER , INTENT(in ) :: jl ! category number 523 524 INTEGER :: ji, jk ! loop indices 525 INTEGER :: ii, ij 526 INTEGER :: numce ! number of points for which conservation is violated 527 REAL(wp) :: meance ! mean conservation error 528 REAL(wp) :: max_cons_err, max_surf_err 529 !!--------------------------------------------------------------------- 530 531 max_cons_err = 1.0_wp ! maximum tolerated conservation error 532 max_surf_err = 0.001_wp ! maximum tolerated surface error 533 534 !-------------------------- 535 ! Increment of energy 536 !-------------------------- 537 ! global 538 DO ji = kideb, kiut 539 dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) + qt_s_fin(ji,jl) - qt_s_in(ji,jl) 540 END DO 541 ! layer by layer 542 dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:) 543 544 !---------------------------------------- 545 ! Atmospheric heat flux, ice heat budget 546 !---------------------------------------- 547 DO ji = kideb, kiut 548 ii = MOD( npb(ji) - 1 , jpi ) + 1 549 ij = ( npb(ji) - 1 ) / jpi + 1 550 fatm (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) 551 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(ii,ij,jl) 552 END DO 553 554 !-------------------- 555 ! Conservation error 556 !-------------------- 557 DO ji = kideb, kiut 558 cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 559 END DO 560 561 numce = 0 562 meance = 0._wp 563 DO ji = kideb, kiut 564 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 565 numce = numce + 1 566 meance = meance + cons_error(ji,jl) 567 ENDIF 568 END DO 569 IF( numce > 0 ) meance = meance / numce 570 571 WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err 572 WRITE(numout,*) ' After lim_thd_dif, category : ', jl 573 WRITE(numout,*) ' Mean conservation error on big error points ', meance, numit 574 WRITE(numout,*) ' Number of points where there is a cons err gt than c.e. : ', numce, numit 575 576 !------------------------------------------------------- 577 ! Surface error due to imbalance between Fatm and Fcsu 578 !------------------------------------------------------- 579 numce = 0 580 meance = 0._wp 581 582 DO ji = kideb, kiut 583 surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) ) 584 IF( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) THEN 585 numce = numce + 1 586 meance = meance + surf_error(ji,jl) 587 ENDIF 588 ENDDO 589 IF( numce > 0 ) meance = meance / numce 590 591 WRITE(numout,*) ' Maximum tolerated surface error : ', max_surf_err 592 WRITE(numout,*) ' After lim_thd_dif, category : ', jl 593 WRITE(numout,*) ' Mean surface error on big error points ', meance, numit 594 WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit 595 596 WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d) 597 WRITE(numout,*) ' fatm : ', fatm(jiindex_1d,jl) 598 WRITE(numout,*) ' t_su : ', t_su_b(jiindex_1d) 599 600 !--------------------------------------- 601 ! Write ice state in case of big errors 602 !--------------------------------------- 603 DO ji = kideb, kiut 604 IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. & 605 ( cons_error(ji,jl) .GT. max_cons_err ) ) THEN 606 ii = MOD( npb(ji) - 1, jpi ) + 1 607 ij = ( npb(ji) - 1 ) / jpi + 1 608 ! 609 WRITE(numout,*) ' alerte 1 ' 610 WRITE(numout,*) ' Untolerated conservation / surface error after ' 611 WRITE(numout,*) ' heat diffusion in the ice ' 612 WRITE(numout,*) ' Category : ', jl 613 WRITE(numout,*) ' ii , ij : ', ii, ij 614 WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij) 615 WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 616 WRITE(numout,*) ' surf_error : ', surf_error(ji,jl) 617 WRITE(numout,*) ' dq_i : ', - dq_i(ji,jl) * r1_rdtice 618 WRITE(numout,*) ' Fdt : ', sum_fluxq(ji,jl) 619 WRITE(numout,*) 620 ! WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl) 621 ! WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl) 622 ! WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) 623 ! WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) 624 ! WRITE(numout,*) ' qt : ', qt_i_fin(ji,jl) + qt_s_fin(ji,jl) 625 WRITE(numout,*) ' ht_i : ', ht_i_b(ji) 626 WRITE(numout,*) ' ht_s : ', ht_s_b(ji) 627 WRITE(numout,*) ' t_su : ', t_su_b(ji) 628 WRITE(numout,*) ' t_s : ', t_s_b(ji,1) 629 WRITE(numout,*) ' t_i : ', t_i_b(ji,1:nlay_i) 630 WRITE(numout,*) ' t_bo : ', t_bo_b(ji) 631 WRITE(numout,*) ' q_i : ', q_i_b(ji,1:nlay_i) 632 WRITE(numout,*) ' s_i : ', s_i_b(ji,1:nlay_i) 633 WRITE(numout,*) ' tmelts : ', rtt - tmut*s_i_b(ji,1:nlay_i) 634 WRITE(numout,*) 635 WRITE(numout,*) ' Fluxes ' 636 WRITE(numout,*) ' ~~~~~~ ' 637 WRITE(numout,*) ' fatm : ', fatm(ji,jl) 638 WRITE(numout,*) ' fc_su : ', fc_su (ji) 639 WRITE(numout,*) ' fstr_inice : ', qsr_ice_1d(ji)*i0(ji) 640 WRITE(numout,*) ' fc_bo : ', - fc_bo_i (ji) 641 WRITE(numout,*) ' foc : ', fbif_1d(ji) 642 WRITE(numout,*) ' fstroc : ', fstroc (ii,ij,jl) 643 WRITE(numout,*) ' i0 : ', i0(ji) 644 WRITE(numout,*) ' qsr_ice : ', (1.0-i0(ji))*qsr_ice_1d(ji) 645 WRITE(numout,*) ' qns_ice : ', qnsr_ice_1d(ji) 646 WRITE(numout,*) ' Conduction fluxes : ' 647 WRITE(numout,*) ' fc_s : ', fc_s(ji,0:nlay_s) 648 WRITE(numout,*) ' fc_i : ', fc_i(ji,0:nlay_i) 649 WRITE(numout,*) 650 WRITE(numout,*) ' Layer by layer ... ' 651 WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) * r1_rdtice 652 WRITE(numout,*) ' dfc_snow : ', fc_s(ji,1) - fc_s(ji,0) 653 DO jk = 1, nlay_i 654 WRITE(numout,*) ' layer : ', jk 655 WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) * r1_rdtice 656 WRITE(numout,*) ' radab : ', radab(ji,jk) 657 WRITE(numout,*) ' dfc_i : ', fc_i(ji,jk) - fc_i(ji,jk-1) 658 WRITE(numout,*) ' tot f : ', fc_i(ji,jk) - fc_i(ji,jk-1) - radab(ji,jk) 659 END DO 660 661 ENDIF 662 ! 663 END DO 664 ! 665 END SUBROUTINE lim_thd_con_dif 666 667 668 SUBROUTINE lim_thd_con_dh( kideb, kiut, jl ) 669 !!----------------------------------------------------------------------- 670 !! *** ROUTINE lim_thd_con_dh *** 671 !! 672 !! ** Purpose : Test energy conservation after enthalpy redistr. 673 !!----------------------------------------------------------------------- 674 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 675 INTEGER, INTENT(in) :: jl ! category number 676 ! 677 INTEGER :: ji ! loop indices 678 INTEGER :: ii, ij, numce ! local integers 679 REAL(wp) :: meance, max_cons_err !local scalar 680 !!--------------------------------------------------------------------- 681 682 max_cons_err = 1._wp 683 684 !-------------------------- 685 ! Increment of energy 686 !-------------------------- 687 DO ji = kideb, kiut 688 dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) + qt_s_fin(ji,jl) - qt_s_in(ji,jl) ! global 689 END DO 690 dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:) ! layer by layer 691 692 !---------------------------------------- 693 ! Atmospheric heat flux, ice heat budget 694 !---------------------------------------- 695 DO ji = kideb, kiut 696 ii = MOD( npb(ji) - 1 , jpi ) + 1 697 ij = ( npb(ji) - 1 ) / jpi + 1 698 699 fatm (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji) ! total heat flux 700 sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(ii,ij,jl) 701 cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 702 END DO 703 704 !-------------------- 705 ! Conservation error 706 !-------------------- 707 DO ji = kideb, kiut 708 cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 709 END DO 710 711 numce = 0 712 meance = 0._wp 713 DO ji = kideb, kiut 714 IF( cons_error(ji,jl) .GT. max_cons_err ) THEN 715 numce = numce + 1 716 meance = meance + cons_error(ji,jl) 717 ENDIF 718 ENDDO 719 IF(numce > 0 ) meance = meance / numce 720 721 WRITE(numout,*) ' Error report - Category : ', jl 722 WRITE(numout,*) ' ~~~~~~~~~~~~ ' 723 WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err 724 WRITE(numout,*) ' After lim_thd_ent, category : ', jl 725 WRITE(numout,*) ' Mean conservation error on big error points ', meance, numit 726 WRITE(numout,*) ' Number of points where there is a cons err gt than 0.1 W/m2 : ', numce, numit 727 728 !--------------------------------------- 729 ! Write ice state in case of big errors 730 !--------------------------------------- 731 DO ji = kideb, kiut 732 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 733 ii = MOD( npb(ji) - 1, jpi ) + 1 734 ij = ( npb(ji) - 1 ) / jpi + 1 735 ! 736 WRITE(numout,*) ' alerte 1 - category : ', jl 737 WRITE(numout,*) ' Untolerated conservation error after limthd_ent ' 738 WRITE(numout,*) ' ii , ij : ', ii, ij 739 WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij) 740 WRITE(numout,*) ' * ' 741 WRITE(numout,*) ' Ftotal : ', sum_fluxq(ji,jl) 742 WRITE(numout,*) ' dq_t : ', - dq_i(ji,jl) * r1_rdtice 743 WRITE(numout,*) ' dq_i : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) * r1_rdtice 744 WRITE(numout,*) ' dq_s : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) * r1_rdtice 745 WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 746 WRITE(numout,*) ' * ' 747 WRITE(numout,*) ' Fluxes --- : ' 748 WRITE(numout,*) ' fatm : ', fatm(ji,jl) 749 WRITE(numout,*) ' foce : ', fbif_1d(ji) 750 WRITE(numout,*) ' fres : ', ftotal_fin(ji) 751 WRITE(numout,*) ' fhbri : ', fhbricat(ii,ij,jl) 752 WRITE(numout,*) ' * ' 753 WRITE(numout,*) ' Heat contents --- : ' 754 WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl) * r1_rdtice 755 WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl) * r1_rdtice 756 WRITE(numout,*) ' qt_in : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) * r1_rdtice 757 WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) * r1_rdtice 758 WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) * r1_rdtice 759 WRITE(numout,*) ' qt_fin : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) * r1_rdtice 760 WRITE(numout,*) ' * ' 761 WRITE(numout,*) ' Ice variables --- : ' 762 WRITE(numout,*) ' ht_i : ', ht_i_b(ji) 763 WRITE(numout,*) ' ht_s : ', ht_s_b(ji) 764 WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji) 765 WRITE(numout,*) ' dh_snowice: ', dh_snowice(ji) 766 WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji) 767 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 768 ENDIF 769 ! 770 END DO 771 ! 772 END SUBROUTINE lim_thd_con_dh 773 774 775 SUBROUTINE lim_thd_enmelt( kideb, kiut ) 776 !!----------------------------------------------------------------------- 777 !! *** ROUTINE lim_thd_enmelt *** 778 !! 779 !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) 421 !! ** Purpose : Computes sea ice temperature (Kelvin) from enthalpy 780 422 !! 781 423 !! ** Method : Formula (Bitz and Lipscomb, 1999) … … 784 426 !! 785 427 INTEGER :: ji, jk ! dummy loop indices 786 REAL(wp) :: ztmelts ! local scalar428 REAL(wp) :: ztmelts, zaaa, zbbb, zccc, zdiscrim ! local scalar 787 429 !!------------------------------------------------------------------- 788 ! 789 DO jk = 1, nlay_i ! Sea ice energy of melting430 ! Recover ice temperature 431 DO jk = 1, nlay_i 790 432 DO ji = kideb, kiut 791 ztmelts = - tmut * s_i_b(ji,jk) + rtt 792 q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) & 793 & + lfus * ( 1.0 - (ztmelts-rtt) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) ) & 794 & - rcp * ( ztmelts-rtt ) ) 795 END DO 433 ztmelts = -tmut * s_i_1d(ji,jk) + rt0 434 ! Conversion q(S,T) -> T (second order equation) 435 zaaa = cpic 436 zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + q_i_1d(ji,jk) * r1_rhoic - lfus 437 zccc = lfus * ( ztmelts - rt0 ) 438 zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 439 t_i_1d(ji,jk) = rt0 - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 440 441 ! mask temperature 442 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 443 t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0 444 END DO 445 END DO 446 447 END SUBROUTINE lim_thd_temp 448 449 SUBROUTINE lim_thd_lam( kideb, kiut ) 450 !!----------------------------------------------------------------------- 451 !! *** ROUTINE lim_thd_lam *** 452 !! 453 !! ** Purpose : Lateral melting in case monocategory 454 !! ( dA = A/2h dh ) 455 !!----------------------------------------------------------------------- 456 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 457 INTEGER :: ji ! dummy loop indices 458 REAL(wp) :: zhi_bef ! ice thickness before thermo 459 REAL(wp) :: zdh_mel, zda_mel ! net melting 460 REAL(wp) :: zvi, zvs ! ice/snow volumes 461 462 DO ji = kideb, kiut 463 zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 464 IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp ) THEN 465 zvi = a_i_1d(ji) * ht_i_1d(ji) 466 zvs = a_i_1d(ji) * ht_s_1d(ji) 467 ! lateral melting = concentration change 468 zhi_bef = ht_i_1d(ji) - zdh_mel 469 rswitch = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) ) 470 zda_mel = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) 471 a_i_1d(ji) = MAX( epsi20, a_i_1d(ji) + zda_mel ) 472 ! adjust thickness 473 ht_i_1d(ji) = zvi / a_i_1d(ji) 474 ht_s_1d(ji) = zvs / a_i_1d(ji) 475 ! retrieve total concentration 476 at_i_1d(ji) = a_i_1d(ji) 477 END IF 796 478 END DO 797 DO jk = 1, nlay_s ! Snow energy of melting 798 DO ji = kideb, kiut 799 q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 800 END DO 801 END DO 802 ! 803 END SUBROUTINE lim_thd_enmelt 479 480 END SUBROUTINE lim_thd_lam 481 482 SUBROUTINE lim_thd_1d2d( nbpb, jl, kn ) 483 !!----------------------------------------------------------------------- 484 !! *** ROUTINE lim_thd_1d2d *** 485 !! 486 !! ** Purpose : move arrays from 1d to 2d and the reverse 487 !!----------------------------------------------------------------------- 488 INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D 489 ! 2= from 1D to 2D 490 INTEGER, INTENT(in) :: nbpb ! size of 1D arrays 491 INTEGER, INTENT(in) :: jl ! ice cat 492 INTEGER :: jk ! dummy loop indices 493 494 SELECT CASE( kn ) 495 496 CASE( 1 ) 497 498 CALL tab_2d_1d( nbpb, at_i_1d (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) ) 499 CALL tab_2d_1d( nbpb, a_i_1d (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 500 CALL tab_2d_1d( nbpb, ht_i_1d (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 501 CALL tab_2d_1d( nbpb, ht_s_1d (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 502 503 CALL tab_2d_1d( nbpb, t_su_1d (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 504 CALL tab_2d_1d( nbpb, sm_i_1d (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 505 DO jk = 1, nlay_s 506 CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 507 CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 508 END DO 509 DO jk = 1, nlay_i 510 CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 511 CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 512 CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 513 END DO 514 515 CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) ) 516 CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 517 CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) ) 518 CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb), fr2_i0 , jpi, jpj, npb(1:nbpb) ) 519 CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 520 CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 521 CALL tab_2d_1d( nbpb, evap_ice_1d (1:nbpb), evap_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 522 CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 523 CALL tab_2d_1d( nbpb, t_bo_1d (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) ) 524 CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) ) 525 CALL tab_2d_1d( nbpb, fhtur_1d (1:nbpb), fhtur , jpi, jpj, npb(1:nbpb) ) 526 CALL tab_2d_1d( nbpb, qlead_1d (1:nbpb), qlead , jpi, jpj, npb(1:nbpb) ) 527 CALL tab_2d_1d( nbpb, fhld_1d (1:nbpb), fhld , jpi, jpj, npb(1:nbpb) ) 528 529 CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw , jpi, jpj, npb(1:nbpb) ) 530 CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub , jpi, jpj, npb(1:nbpb) ) 531 532 CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog , jpi, jpj, npb(1:nbpb) ) 533 CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom , jpi, jpj, npb(1:nbpb) ) 534 CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum , jpi, jpj, npb(1:nbpb) ) 535 CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni , jpi, jpj, npb(1:nbpb) ) 536 CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res , jpi, jpj, npb(1:nbpb) ) 537 CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr , jpi, jpj, npb(1:nbpb) ) 538 539 CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog , jpi, jpj, npb(1:nbpb) ) 540 CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom , jpi, jpj, npb(1:nbpb) ) 541 CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum , jpi, jpj, npb(1:nbpb) ) 542 CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni , jpi, jpj, npb(1:nbpb) ) 543 CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) ) 544 CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res , jpi, jpj, npb(1:nbpb) ) 545 546 CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd , jpi, jpj, npb(1:nbpb) ) 547 CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr , jpi, jpj, npb(1:nbpb) ) 548 CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum , jpi, jpj, npb(1:nbpb) ) 549 CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom , jpi, jpj, npb(1:nbpb) ) 550 CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog , jpi, jpj, npb(1:nbpb) ) 551 CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif , jpi, jpj, npb(1:nbpb) ) 552 CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw , jpi, jpj, npb(1:nbpb) ) 553 CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw , jpi, jpj, npb(1:nbpb) ) 554 CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub , jpi, jpj, npb(1:nbpb) ) 555 CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err , jpi, jpj, npb(1:nbpb) ) 556 CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res , jpi, jpj, npb(1:nbpb) ) 557 CALL tab_2d_1d( nbpb, hfx_err_dif_1d (1:nbpb), hfx_err_dif , jpi, jpj, npb(1:nbpb) ) 558 CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) ) 559 560 CASE( 2 ) 561 562 CALL tab_1d_2d( nbpb, at_i , npb, at_i_1d (1:nbpb) , jpi, jpj ) 563 CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_1d (1:nbpb) , jpi, jpj ) 564 CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_1d (1:nbpb) , jpi, jpj ) 565 CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_1d (1:nbpb) , jpi, jpj ) 566 CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_1d (1:nbpb) , jpi, jpj ) 567 CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_1d (1:nbpb) , jpi, jpj ) 568 DO jk = 1, nlay_s 569 CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d (1:nbpb,jk), jpi, jpj) 570 CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d (1:nbpb,jk), jpi, jpj) 571 END DO 572 DO jk = 1, nlay_i 573 CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d (1:nbpb,jk), jpi, jpj) 574 CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d (1:nbpb,jk), jpi, jpj) 575 CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d (1:nbpb,jk), jpi, jpj) 576 END DO 577 CALL tab_1d_2d( nbpb, qlead , npb, qlead_1d (1:nbpb) , jpi, jpj ) 578 579 CALL tab_1d_2d( nbpb, wfx_snw , npb, wfx_snw_1d(1:nbpb) , jpi, jpj ) 580 CALL tab_1d_2d( nbpb, wfx_sub , npb, wfx_sub_1d(1:nbpb) , jpi, jpj ) 581 582 CALL tab_1d_2d( nbpb, wfx_bog , npb, wfx_bog_1d(1:nbpb) , jpi, jpj ) 583 CALL tab_1d_2d( nbpb, wfx_bom , npb, wfx_bom_1d(1:nbpb) , jpi, jpj ) 584 CALL tab_1d_2d( nbpb, wfx_sum , npb, wfx_sum_1d(1:nbpb) , jpi, jpj ) 585 CALL tab_1d_2d( nbpb, wfx_sni , npb, wfx_sni_1d(1:nbpb) , jpi, jpj ) 586 CALL tab_1d_2d( nbpb, wfx_res , npb, wfx_res_1d(1:nbpb) , jpi, jpj ) 587 CALL tab_1d_2d( nbpb, wfx_spr , npb, wfx_spr_1d(1:nbpb) , jpi, jpj ) 588 589 CALL tab_1d_2d( nbpb, sfx_bog , npb, sfx_bog_1d(1:nbpb) , jpi, jpj ) 590 CALL tab_1d_2d( nbpb, sfx_bom , npb, sfx_bom_1d(1:nbpb) , jpi, jpj ) 591 CALL tab_1d_2d( nbpb, sfx_sum , npb, sfx_sum_1d(1:nbpb) , jpi, jpj ) 592 CALL tab_1d_2d( nbpb, sfx_sni , npb, sfx_sni_1d(1:nbpb) , jpi, jpj ) 593 CALL tab_1d_2d( nbpb, sfx_res , npb, sfx_res_1d(1:nbpb) , jpi, jpj ) 594 CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj ) 595 596 CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj ) 597 CALL tab_1d_2d( nbpb, hfx_spr , npb, hfx_spr_1d(1:nbpb) , jpi, jpj ) 598 CALL tab_1d_2d( nbpb, hfx_sum , npb, hfx_sum_1d(1:nbpb) , jpi, jpj ) 599 CALL tab_1d_2d( nbpb, hfx_bom , npb, hfx_bom_1d(1:nbpb) , jpi, jpj ) 600 CALL tab_1d_2d( nbpb, hfx_bog , npb, hfx_bog_1d(1:nbpb) , jpi, jpj ) 601 CALL tab_1d_2d( nbpb, hfx_dif , npb, hfx_dif_1d(1:nbpb) , jpi, jpj ) 602 CALL tab_1d_2d( nbpb, hfx_opw , npb, hfx_opw_1d(1:nbpb) , jpi, jpj ) 603 CALL tab_1d_2d( nbpb, hfx_snw , npb, hfx_snw_1d(1:nbpb) , jpi, jpj ) 604 CALL tab_1d_2d( nbpb, hfx_sub , npb, hfx_sub_1d(1:nbpb) , jpi, jpj ) 605 CALL tab_1d_2d( nbpb, hfx_err , npb, hfx_err_1d(1:nbpb) , jpi, jpj ) 606 CALL tab_1d_2d( nbpb, hfx_res , npb, hfx_res_1d(1:nbpb) , jpi, jpj ) 607 CALL tab_1d_2d( nbpb, hfx_err_rem , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj ) 608 CALL tab_1d_2d( nbpb, hfx_err_dif , npb, hfx_err_dif_1d(1:nbpb), jpi, jpj ) 609 ! 610 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 611 CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 612 ! 613 END SELECT 614 615 END SUBROUTINE lim_thd_1d2d 804 616 805 617 … … 817 629 !!------------------------------------------------------------------- 818 630 INTEGER :: ios ! Local integer output status for namelist read 819 NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, & 820 & hicmin, hiclim, & 821 & sbeta , parlat, hakspl, hibspl, exld, & 822 & hakdif, hnzst , thth , parsub, alphs, betas, & 823 & kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi 631 NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, & 632 & rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 633 & nn_monocat, ln_it_qnsice 824 634 !!------------------------------------------------------------------- 825 635 ! … … 839 649 IF(lwm) WRITE ( numoni, namicethd ) 840 650 ! 651 IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN 652 nn_monocat = 0 653 IF(lwp) WRITE(numout, *) ' nn_monocat must be 0 in multi-category case ' 654 ENDIF 655 656 ! 841 657 IF(lwp) THEN ! control print 842 658 WRITE(numout,*) 843 659 WRITE(numout,*)' Namelist of ice parameters for ice thermodynamic computation ' 844 WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt 845 WRITE(numout,*)' ice thick. for lateral accretion in NH (SH) hiccrit(1/2) = ', hiccrit 846 WRITE(numout,*)' Frazil ice thickness as a function of wind or not fraz_swi = ', fraz_swi 847 WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom maxfrazb = ', maxfrazb 848 WRITE(numout,*)' Thresold relative drift speed for collection of frazil vfrazb = ', vfrazb 849 WRITE(numout,*)' Squeezing coefficient for collection of frazil Cfrazb = ', Cfrazb 850 WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin 851 WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim 660 WRITE(numout,*)' ice thick. for lateral accretion rn_hnewice = ', rn_hnewice 661 WRITE(numout,*)' Frazil ice thickness as a function of wind or not ln_frazil = ', ln_frazil 662 WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom rn_maxfrazb = ', rn_maxfrazb 663 WRITE(numout,*)' Thresold relative drift speed for collection of frazil rn_vfrazb = ', rn_vfrazb 664 WRITE(numout,*)' Squeezing coefficient for collection of frazil rn_Cfrazb = ', rn_Cfrazb 665 WRITE(numout,*)' minimum ice thickness rn_himin = ', rn_himin 852 666 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 853 WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta 854 WRITE(numout,*)' percentage of energy used for lateral ablation parlat = ', parlat 855 WRITE(numout,*)' slope of distr. for Hakkinen-Mellor lateral melting hakspl = ', hakspl 856 WRITE(numout,*)' slope of distribution for Hibler lateral melting hibspl = ', hibspl 857 WRITE(numout,*)' exponent for leads-closure rate exld = ', exld 858 WRITE(numout,*)' coefficient for diffusions of ice and snow hakdif = ', hakdif 859 WRITE(numout,*)' threshold thick. for comp. of eq. thermal conductivity zhth = ', thth 860 WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst 861 WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub 862 WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs 863 WRITE(numout,*)' coefficient for ice-lead partition of snowfall betas = ', betas 864 WRITE(numout,*)' extinction radiation parameter in sea ice (1.0) kappa_i = ', kappa_i 865 WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nconv_i_thd = ', nconv_i_thd 866 WRITE(numout,*)' maximal err. on T for heat diffusion computation maxer_i_thd = ', maxer_i_thd 867 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice thcon_i_swi = ', thcon_i_swi 667 WRITE(numout,*)' coefficient for ice-lead partition of snowfall rn_betas = ', rn_betas 668 WRITE(numout,*)' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i 669 WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nn_conv_dif = ', nn_conv_dif 670 WRITE(numout,*)' maximal err. on T for heat diffusion computation rn_terr_dif = ', rn_terr_dif 671 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice nn_ice_thcon = ', nn_ice_thcon 672 WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i 673 WRITE(numout,*)' virtual ITD mono-category parameterizations (1) or not nn_monocat = ', nn_monocat 674 WRITE(numout,*)' iterate the surface non-solar flux (T) or not (F) ln_it_qnsice = ', ln_it_qnsice 868 675 ENDIF 869 !870 rcdsn = hakdif * rcdsn871 rcdic = hakdif * rcdic872 676 ! 873 677 END SUBROUTINE lim_thd_init
Note: See TracChangeset
for help on using the changeset viewer.
