- Timestamp:
- 2015-02-11T11:50:34+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4624 r5075 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, oatte24 USE oce , ONLY : fraqsr_1lev 25 25 USE ice ! LIM: sea-ice variables 26 26 USE par_ice ! LIM: sea-ice parameters … … 43 43 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 44 44 USE timing ! Timing 45 USE limcons ! conservation tests 45 46 46 47 IMPLICIT NONE … … 49 50 PUBLIC lim_thd ! called by limstp module 50 51 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 !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 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 thermo. cal. 86 INTEGER :: ii, ij ! temporary dummy loop index 87 REAL(wp) :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 88 REAL(wp) :: zch = 0.0057_wp ! heat transfer coefficient 89 REAL(wp) :: zareamin 90 REAL(wp) :: zfric_u, zqld, zqfr 91 ! 92 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 93 ! 94 REAL(wp), POINTER, DIMENSION(:,:) :: zqsr, zqns 95 95 !!------------------------------------------------------------------- 96 CALL wrk_alloc( jpi, jpj, zqsr, zqns ) 97 96 98 IF( nn_timing == 1 ) CALL timing_start('limthd') 97 99 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 !------------------------------------------------------------------------------! 100 ! conservation test 101 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 102 103 !------------------------------------------------------------------------! 104 ! 1) Initialization of some variables ! 105 !------------------------------------------------------------------------! 106 ftr_ice(:,:,:) = 0._wp ! part of solar radiation transmitted through the ice 107 114 108 115 109 !-------------------- … … 121 115 DO jj = 1, jpj 122 116 DO ji = 1, jpi 117 !0 if no ice and 1 if yes 118 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) 123 119 !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 120 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 121 !convert units ! very important that this line is here 122 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 129 123 END DO 130 124 END DO … … 133 127 DO jj = 1, jpj 134 128 DO ji = 1, jpi 129 !0 if no ice and 1 if yes 130 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 ) ) 135 131 !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 ) ) 132 e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s ) 139 133 !convert units ! very important that this line is here 140 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb134 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac 141 135 END DO 142 136 END DO 143 137 END DO 144 138 END DO 145 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 139 156 140 ! 2) Partial computation of forcing for the thermodynamic sea ice model. ! 157 141 !-----------------------------------------------------------------------------! 142 143 !--- Ocean solar and non solar fluxes to be used in zqld 144 IF ( .NOT. lk_cpl ) THEN ! --- forced case, fluxes to the lead are the same as over the ocean 145 ! 146 zqsr(:,:) = qsr(:,:) ; zqns(:,:) = qns(:,:) 147 ! 148 ELSE ! --- coupled case, fluxes to the lead are total - intercepted 149 ! 150 zqsr(:,:) = qsr_tot(:,:) ; zqns(:,:) = qns_tot(:,:) 151 ! 152 DO jl = 1, jpl 153 DO jj = 1, jpj 154 DO ji = 1, jpi 155 zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 156 zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 157 END DO 158 END DO 159 END DO 160 ! 161 ENDIF 158 162 159 163 !CDIR NOVERRCHK … … 161 165 !CDIR NOVERRCHK 162 166 DO ji = 1, jpi 163 zinda = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) + epsi10 ) ) )167 rswitch = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 164 168 ! 165 169 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 168 172 ! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean 169 173 ! ! 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 174 ! 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 ) ) 175 176 ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- ! 177 ! REMARK valid at least in forced mode from clem 178 ! precip is included in qns but not in qns_ice 179 IF ( lk_cpl ) THEN 180 zqld = tms(ji,jj) * rdt_ice * & 181 & ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) & ! pfrld already included in coupled mode 182 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip 183 & ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 184 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) 185 ELSE 186 zqld = tms(ji,jj) * rdt_ice * & 187 & ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) ) & 188 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip 189 & ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 190 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) 191 ENDIF 192 193 !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 194 zqfr = tms(ji,jj) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 195 196 !-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice 197 qlead(ji,jj) = MIN( 0._wp , zqld - zqfr ) 198 199 ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting 200 IF( at_i(ji,jj) > epsi10 .AND. zqld > 0._wp ) THEN 201 fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90 202 qlead(ji,jj) = 0._wp 203 ELSE 204 fhld (ji,jj) = 0._wp 205 ENDIF 194 206 ! 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 ! 207 !-- Energy from the turbulent oceanic heat flux --- ! 208 !clem zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) 209 zfric_u = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 210 fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2 211 ! upper bound for fhtur: we do not want SST to drop below Tfreeze. 212 ! So we say that the heat retrieved from the ocean (fhtur+fhld) must be < to the heat necessary to reach Tfreeze (zqfr) 213 ! This is not a clean budget, so that should be corrected at some point 214 fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - fhld(ji,jj) - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 215 216 ! ----------------------------------------- 217 ! Net heat flux on top of ice-ocean [W.m-2] 218 ! ----------------------------------------- 219 ! First step here : heat flux at the ocean surface + precip 220 ! Second step below : heat flux at the ice surface (after limthd_dif) 221 hfx_in(ji,jj) = hfx_in(ji,jj) & 222 ! heat flux above the ocean 223 & + pfrld(ji,jj) * ( zqns(ji,jj) + zqsr(ji,jj) ) & 224 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 225 & + ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 226 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) 227 228 ! ----------------------------------------------------------------------------- 229 ! Net heat flux that is retroceded to the ocean or taken from the ocean [W.m-2] 230 ! ----------------------------------------------------------------------------- 231 ! First step here : non solar + precip - qlead - qturb 232 ! Second step in limthd_dh : heat remaining if total melt (zq_rema) 233 ! Third step in limsbc : heat from ice-ocean mass exchange (zf_mass) + solar 234 hfx_out(ji,jj) = hfx_out(ji,jj) & 235 ! Non solar heat flux received by the ocean 236 & + pfrld(ji,jj) * qns(ji,jj) & 237 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 238 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) & 239 & * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 240 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) & 241 ! heat flux taken from the ocean where there is open water ice formation 242 & - qlead(ji,jj) * r1_rdtice & 243 ! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth) 244 & - at_i(ji,jj) * fhtur(ji,jj) & 245 & - at_i(ji,jj) * fhld(ji,jj) 246 205 247 END DO 206 248 END DO … … 234 276 DO jj = mj0(jjindx), mj1(jjindx) 235 277 jiindex_1d = (jj - 1) * jpi + ji 278 WRITE(numout,*) ' lim_thd : Category no : ', jl 236 279 END DO 237 280 END DO … … 250 293 !------------------------- 251 294 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) )295 CALL tab_2d_1d( nbpb, at_i_1d (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) ) 296 CALL tab_2d_1d( nbpb, a_i_1d (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 297 CALL tab_2d_1d( nbpb, ht_i_1d (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 298 CALL tab_2d_1d( nbpb, ht_s_1d (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 299 300 CALL tab_2d_1d( nbpb, t_su_1d (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 301 CALL tab_2d_1d( nbpb, sm_i_1d (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 259 302 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) )303 CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 304 CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 262 305 END DO 263 306 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) )307 CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 308 CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 309 CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 267 310 END DO 268 311 … … 271 314 CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) ) 272 315 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 316 CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 317 CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 318 IF( .NOT. lk_cpl ) THEN 319 CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 320 CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 321 ENDIF 278 322 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) )323 CALL tab_2d_1d( nbpb, t_bo_1d (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) ) 280 324 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) ) 325 CALL tab_2d_1d( nbpb, fhtur_1d (1:nbpb), fhtur , jpi, jpj, npb(1:nbpb) ) 326 CALL tab_2d_1d( nbpb, qlead_1d (1:nbpb), qlead , jpi, jpj, npb(1:nbpb) ) 327 CALL tab_2d_1d( nbpb, fhld_1d (1:nbpb), fhld , jpi, jpj, npb(1:nbpb) ) 328 329 CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw , jpi, jpj, npb(1:nbpb) ) 330 CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub , jpi, jpj, npb(1:nbpb) ) 331 332 CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog , jpi, jpj, npb(1:nbpb) ) 333 CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom , jpi, jpj, npb(1:nbpb) ) 334 CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum , jpi, jpj, npb(1:nbpb) ) 335 CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni , jpi, jpj, npb(1:nbpb) ) 336 CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res , jpi, jpj, npb(1:nbpb) ) 337 CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr , jpi, jpj, npb(1:nbpb) ) 338 339 CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog , jpi, jpj, npb(1:nbpb) ) 340 CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom , jpi, jpj, npb(1:nbpb) ) 341 CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum , jpi, jpj, npb(1:nbpb) ) 342 CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni , jpi, jpj, npb(1:nbpb) ) 289 343 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 344 CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res , jpi, jpj, npb(1:nbpb) ) 345 346 CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd , jpi, jpj, npb(1:nbpb) ) 347 CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr , jpi, jpj, npb(1:nbpb) ) 348 CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum , jpi, jpj, npb(1:nbpb) ) 349 CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom , jpi, jpj, npb(1:nbpb) ) 350 CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog , jpi, jpj, npb(1:nbpb) ) 351 CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif , jpi, jpj, npb(1:nbpb) ) 352 CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw , jpi, jpj, npb(1:nbpb) ) 353 CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw , jpi, jpj, npb(1:nbpb) ) 354 CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub , jpi, jpj, npb(1:nbpb) ) 355 CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err , jpi, jpj, npb(1:nbpb) ) 356 CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res , jpi, jpj, npb(1:nbpb) ) 357 CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) ) 358 296 359 !-------------------------------- 297 360 ! 4.3) Thermodynamic processes 298 361 !-------------------------------- 299 362 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 ) 363 !---------------------------------! 364 ! Ice/Snow Temperature profile ! 365 !---------------------------------! 366 CALL lim_thd_dif( 1, nbpb ) 367 368 !---------------------------------! 369 ! Ice/Snow thicnkess ! 370 !---------------------------------! 371 CALL lim_thd_dh( 1, nbpb ) 372 373 ! --- Ice enthalpy remapping --- ! 374 CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) 375 376 !---------------------------------! 377 ! --- Ice salinity --- ! 378 !---------------------------------! 379 CALL lim_thd_sal( 1, nbpb ) 380 381 !---------------------------------! 382 ! --- temperature update --- ! 383 !---------------------------------! 384 CALL lim_thd_temp( 1, nbpb ) 327 385 328 386 !-------------------------------- … … 330 388 !-------------------------------- 331 389 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 )390 CALL tab_1d_2d( nbpb, at_i , npb, at_i_1d (1:nbpb) , jpi, jpj ) 391 CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_1d (1:nbpb) , jpi, jpj ) 392 CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_1d (1:nbpb) , jpi, jpj ) 393 CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_1d (1:nbpb) , jpi, jpj ) 394 CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_1d (1:nbpb) , jpi, jpj ) 395 CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_1d (1:nbpb) , jpi, jpj ) 338 396 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)397 CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d (1:nbpb,jk), jpi, jpj) 398 CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d (1:nbpb,jk), jpi, jpj) 341 399 END DO 342 400 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 ) 401 CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d (1:nbpb,jk), jpi, jpj) 402 CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d (1:nbpb,jk), jpi, jpj) 403 CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d (1:nbpb,jk), jpi, jpj) 404 END DO 405 CALL tab_1d_2d( nbpb, qlead , npb, qlead_1d (1:nbpb) , jpi, jpj ) 406 407 CALL tab_1d_2d( nbpb, wfx_snw , npb, wfx_snw_1d(1:nbpb) , jpi, jpj ) 408 CALL tab_1d_2d( nbpb, wfx_sub , npb, wfx_sub_1d(1:nbpb) , jpi, jpj ) 409 410 CALL tab_1d_2d( nbpb, wfx_bog , npb, wfx_bog_1d(1:nbpb) , jpi, jpj ) 411 CALL tab_1d_2d( nbpb, wfx_bom , npb, wfx_bom_1d(1:nbpb) , jpi, jpj ) 412 CALL tab_1d_2d( nbpb, wfx_sum , npb, wfx_sum_1d(1:nbpb) , jpi, jpj ) 413 CALL tab_1d_2d( nbpb, wfx_sni , npb, wfx_sni_1d(1:nbpb) , jpi, jpj ) 414 CALL tab_1d_2d( nbpb, wfx_res , npb, wfx_res_1d(1:nbpb) , jpi, jpj ) 415 CALL tab_1d_2d( nbpb, wfx_spr , npb, wfx_spr_1d(1:nbpb) , jpi, jpj ) 416 417 CALL tab_1d_2d( nbpb, sfx_bog , npb, sfx_bog_1d(1:nbpb) , jpi, jpj ) 418 CALL tab_1d_2d( nbpb, sfx_bom , npb, sfx_bom_1d(1:nbpb) , jpi, jpj ) 419 CALL tab_1d_2d( nbpb, sfx_sum , npb, sfx_sum_1d(1:nbpb) , jpi, jpj ) 420 CALL tab_1d_2d( nbpb, sfx_sni , npb, sfx_sni_1d(1:nbpb) , jpi, jpj ) 421 CALL tab_1d_2d( nbpb, sfx_res , npb, sfx_res_1d(1:nbpb) , jpi, jpj ) 422 CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj ) 423 424 CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj ) 425 CALL tab_1d_2d( nbpb, hfx_spr , npb, hfx_spr_1d(1:nbpb) , jpi, jpj ) 426 CALL tab_1d_2d( nbpb, hfx_sum , npb, hfx_sum_1d(1:nbpb) , jpi, jpj ) 427 CALL tab_1d_2d( nbpb, hfx_bom , npb, hfx_bom_1d(1:nbpb) , jpi, jpj ) 428 CALL tab_1d_2d( nbpb, hfx_bog , npb, hfx_bog_1d(1:nbpb) , jpi, jpj ) 429 CALL tab_1d_2d( nbpb, hfx_dif , npb, hfx_dif_1d(1:nbpb) , jpi, jpj ) 430 CALL tab_1d_2d( nbpb, hfx_opw , npb, hfx_opw_1d(1:nbpb) , jpi, jpj ) 431 CALL tab_1d_2d( nbpb, hfx_snw , npb, hfx_snw_1d(1:nbpb) , jpi, jpj ) 432 CALL tab_1d_2d( nbpb, hfx_sub , npb, hfx_sub_1d(1:nbpb) , jpi, jpj ) 433 CALL tab_1d_2d( nbpb, hfx_err , npb, hfx_err_1d(1:nbpb) , jpi, jpj ) 434 CALL tab_1d_2d( nbpb, hfx_res , npb, hfx_res_1d(1:nbpb) , jpi, jpj ) 435 CALL tab_1d_2d( nbpb, hfx_err_rem , npb, hfx_err_rem_1d(1:nbpb) , jpi, jpj ) 358 436 ! 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 !+++++ 437 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 438 CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 373 439 ! 374 440 IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? … … 384 450 ! 5.1) Ice heat content 385 451 !------------------------ 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 ) ) 452 ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 388 453 DO jl = 1, jpl 389 454 DO jk = 1, nlay_i 390 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) * zcoef455 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) / ( unit_fac * REAL( nlay_i ) ) 391 456 END DO 392 457 END DO … … 395 460 ! 5.2) Snow heat content 396 461 !------------------------ 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 ) ) 462 ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 399 463 DO jl = 1, jpl 400 464 DO jk = 1, nlay_s 401 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) * zcoef465 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) / ( unit_fac * REAL( nlay_s ) ) 402 466 END DO 403 467 END DO … … 411 475 ! 5.4) Diagnostic thermodynamic growth rates 412 476 !-------------------------------------------- 413 !clem@useless d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes414 !clem@mv-to-itd dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday415 416 IF( con_i .AND. jiindex_1d > 0 ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)417 418 477 IF(ln_ctl) THEN ! Control print 419 478 CALL prt_ctl_info(' ') … … 448 507 ENDIF 449 508 ! 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 455 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 ) 509 ! 510 CALL wrk_dealloc( jpi, jpj, zqsr, zqns ) 511 512 ! 513 ! conservation test 514 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 475 515 ! 476 516 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)517 518 END SUBROUTINE lim_thd 519 520 SUBROUTINE lim_thd_temp( kideb, kiut ) 481 521 !!----------------------------------------------------------------------- 482 !! *** ROUTINE lim_thd_ glohec***522 !! *** ROUTINE lim_thd_temp *** 483 523 !! 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) 524 !! ** Purpose : Computes sea ice temperature (Kelvin) from enthalpy 780 525 !! 781 526 !! ** Method : Formula (Bitz and Lipscomb, 1999) … … 784 529 !! 785 530 INTEGER :: ji, jk ! dummy loop indices 786 REAL(wp) :: ztmelts ! local scalar531 REAL(wp) :: ztmelts, zaaa, zbbb, zccc, zdiscrim ! local scalar 787 532 !!------------------------------------------------------------------- 788 ! 789 DO jk = 1, nlay_i ! Sea ice energy of melting533 ! Recover ice temperature 534 DO jk = 1, nlay_i 790 535 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 796 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 804 536 ztmelts = -tmut * s_i_1d(ji,jk) + rtt 537 ! Conversion q(S,T) -> T (second order equation) 538 zaaa = cpic 539 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) / rhoic - lfus 540 zccc = lfus * ( ztmelts - rtt ) 541 zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 542 t_i_1d(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 543 544 ! mask temperature 545 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 546 t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rtt 547 END DO 548 END DO 549 550 END SUBROUTINE lim_thd_temp 805 551 806 552 SUBROUTINE lim_thd_init … … 818 564 INTEGER :: ios ! Local integer output status for namelist read 819 565 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, & 566 & hiclim, hnzst, parsub, betas, & 823 567 & kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi 824 568 !!------------------------------------------------------------------- … … 838 582 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp ) 839 583 IF(lwm) WRITE ( numoni, namicethd ) 584 585 IF( lk_cpl .AND. parsub /= 0.0 ) CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' ) 840 586 ! 841 587 IF(lwp) THEN ! control print … … 843 589 WRITE(numout,*)' Namelist of ice parameters for ice thermodynamic computation ' 844 590 WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt 845 WRITE(numout,*)' ice thick. for lateral accretion in NH (SH) hiccrit(1/2)= ', hiccrit591 WRITE(numout,*)' ice thick. for lateral accretion hiccrit = ', hiccrit 846 592 WRITE(numout,*)' Frazil ice thickness as a function of wind or not fraz_swi = ', fraz_swi 847 593 WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom maxfrazb = ', maxfrazb 848 594 WRITE(numout,*)' Thresold relative drift speed for collection of frazil vfrazb = ', vfrazb 849 595 WRITE(numout,*)' Squeezing coefficient for collection of frazil Cfrazb = ', Cfrazb 850 WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin851 596 WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim 852 597 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 853 WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta854 WRITE(numout,*)' percentage of energy used for lateral ablation parlat = ', parlat855 WRITE(numout,*)' slope of distr. for Hakkinen-Mellor lateral melting hakspl = ', hakspl856 WRITE(numout,*)' slope of distribution for Hibler lateral melting hibspl = ', hibspl857 WRITE(numout,*)' exponent for leads-closure rate exld = ', exld858 WRITE(numout,*)' coefficient for diffusions of ice and snow hakdif = ', hakdif859 WRITE(numout,*)' threshold thick. for comp. of eq. thermal conductivity zhth = ', thth860 598 WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst 861 599 WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub 862 WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs863 600 WRITE(numout,*)' coefficient for ice-lead partition of snowfall betas = ', betas 864 601 WRITE(numout,*)' extinction radiation parameter in sea ice (1.0) kappa_i = ', kappa_i … … 866 603 WRITE(numout,*)' maximal err. on T for heat diffusion computation maxer_i_thd = ', maxer_i_thd 867 604 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice thcon_i_swi = ', thcon_i_swi 605 WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i 868 606 ENDIF 869 !870 rcdsn = hakdif * rcdsn871 rcdic = hakdif * rcdic872 607 ! 873 608 END SUBROUTINE lim_thd_init
Note: See TracChangeset
for help on using the changeset viewer.
