Changeset 6515 for branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
- Timestamp:
- 2016-05-09T16:42:28+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r6469 r6515 29 29 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 30 30 USE wrk_nemo ! work arrays 31 USE fldread ! read input fields 32 USE iom 31 33 32 34 IMPLICIT NONE … … 35 37 PUBLIC lim_istate ! routine called by lim_init.F90 36 38 37 ! !!** init namelist (namiceini) ** 38 REAL(wp) :: rn_thres_sst ! threshold water temperature for initial sea ice 39 REAL(wp) :: rn_hts_ini_n ! initial snow thickness in the north 40 REAL(wp) :: rn_hts_ini_s ! initial snow thickness in the south 41 REAL(wp) :: rn_hti_ini_n ! initial ice thickness in the north 42 REAL(wp) :: rn_hti_ini_s ! initial ice thickness in the south 43 REAL(wp) :: rn_ati_ini_n ! initial leads area in the north 44 REAL(wp) :: rn_ati_ini_s ! initial leads area in the south 45 REAL(wp) :: rn_smi_ini_n ! initial salinity 46 REAL(wp) :: rn_smi_ini_s ! initial salinity 47 REAL(wp) :: rn_tmi_ini_n ! initial temperature 48 REAL(wp) :: rn_tmi_ini_s ! initial temperature 49 50 LOGICAL :: ln_iceini ! initialization or not 39 INTEGER , PARAMETER :: jpfldi = 6 ! maximum number of files to read 40 INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m) at T-point 41 INTEGER , PARAMETER :: jp_hts = 2 ! index of snow thicknes (m) at T-point 42 INTEGER , PARAMETER :: jp_ati = 3 ! index of ice fraction (%) at T-point 43 INTEGER , PARAMETER :: jp_tsu = 4 ! index of ice surface temp (K) at T-point 44 INTEGER , PARAMETER :: jp_tmi = 5 ! index of ice temp at T-point 45 INTEGER , PARAMETER :: jp_smi = 6 ! index of ice sali at T-point 46 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read) 51 47 !!---------------------------------------------------------------------- 52 48 !! LIM 3.0, UCL-LOCEAN-IPSL (2008) … … 89 85 REAL(wp) :: ztmelts, zdh 90 86 INTEGER :: i_hemis, i_fill, jl0 91 REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 92 REAL(wp), POINTER, DIMENSION(:) :: zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini 93 REAL(wp), POINTER, DIMENSION(:,:) :: zh_i_ini, za_i_ini, zv_i_ini 87 REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 94 88 REAL(wp), POINTER, DIMENSION(:,:) :: zswitch ! ice indicator 95 INTEGER, POINTER, DIMENSION(:,:) :: zhemis ! hemispheric index 96 !-------------------------------------------------------------------- 97 98 CALL wrk_alloc( jpi, jpj, zswitch ) 99 CALL wrk_alloc( jpi, jpj, zhemis ) 100 CALL wrk_alloc( jpl, 2, zh_i_ini, za_i_ini, zv_i_ini ) 101 CALL wrk_alloc( 2, zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 89 REAL(wp), POINTER, DIMENSION(:,:) :: zht_i_ini, zat_i_ini, zvt_i_ini !data from namelist or nc file 90 REAL(wp), POINTER, DIMENSION(:,:) :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 91 REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini, zv_i_ini !data by cattegories to fill 92 !-------------------------------------------------------------------- 93 94 CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini, za_i_ini, zv_i_ini ) 95 CALL wrk_alloc( jpi, jpj, zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 96 CALL wrk_alloc( jpi, jpj, zswitch ) 102 97 103 98 IF(lwp) WRITE(numout,*) … … 121 116 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 122 117 123 IF( ln_iceini ) THEN 124 125 !-------------------------------------------------------------------- 126 ! 2) Basal temperature, ice mask and hemispheric index 127 !-------------------------------------------------------------------- 128 129 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 130 DO ji = 1, jpi 131 IF( ( sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 132 zswitch(ji,jj) = 0._wp * tmask(ji,jj,1) ! no ice 133 ELSE 134 zswitch(ji,jj) = 1._wp * tmask(ji,jj,1) ! ice 135 ENDIF 118 119 IF( ln_limini ) THEN 120 121 !-------------------------------------------------------------------- 122 ! 2) Basal temperature, ice mask and hemispheric index 123 !-------------------------------------------------------------------- 124 125 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 126 DO ji = 1, jpi 127 IF( ( sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 128 zswitch(ji,jj) = 0._wp * tmask(ji,jj,1) ! no ice 129 ELSE 130 zswitch(ji,jj) = 1._wp * tmask(ji,jj,1) ! ice 131 ENDIF 132 END DO 136 133 END DO 137 END DO 138 139 140 ! Hemispheric index 141 DO jj = 1, jpj 142 DO ji = 1, jpi 143 IF( fcor(ji,jj) >= 0._wp ) THEN 144 zhemis(ji,jj) = 1 ! Northern hemisphere 145 ELSE 146 zhemis(ji,jj) = 2 ! Southern hemisphere 147 ENDIF 148 END DO 149 END DO 150 151 !-------------------------------------------------------------------- 152 ! 3) Initialization of sea ice state variables 153 !-------------------------------------------------------------------- 154 155 !----------------------------- 156 ! 3.1) Hemisphere-dependent arrays 157 !----------------------------- 158 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 159 zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s ! ice thickness 160 zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s ! snow depth 161 zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s ! ice concentration 162 zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s ! bulk ice salinity 163 ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s ! temperature (ice and snow) 164 165 zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:) ! ice volume 166 167 !--------------------------------------------------------------------- 168 ! 3.2) Distribute ice concentration and thickness into the categories 169 !--------------------------------------------------------------------- 170 ! a gaussian distribution for ice concentration is used 171 ! then we check whether the distribution fullfills 172 ! volume and area conservation, positivity and ice categories bounds 173 DO i_hemis = 1, 2 174 175 ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 176 177 ! note for the great nemo engineers: 178 ! only very few of the WRITE statements are necessary for the reference version 179 ! they were one day useful, but now i personally doubt of their 180 ! potential for bringing anything useful 181 182 DO i_fill = jpl, 1, -1 183 IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 184 !---------------------------- 185 ! fill the i_fill categories 186 !---------------------------- 187 ! *** 1 category to fill 188 IF ( i_fill .EQ. 1 ) THEN 189 zh_i_ini(1,i_hemis) = zht_i_ini(i_hemis) 190 za_i_ini(1,i_hemis) = zat_i_ini(i_hemis) 191 zh_i_ini(2:jpl,i_hemis) = 0._wp 192 za_i_ini(2:jpl,i_hemis) = 0._wp 193 ELSE 194 195 ! *** >1 categores to fill 196 !--- Ice thicknesses in the i_fill - 1 first categories 197 DO jl = 1, i_fill - 1 198 zh_i_ini(jl,i_hemis) = hi_mean(jl) 199 END DO 200 201 !--- jl0: most likely index where cc will be maximum 202 DO jl = 1, jpl 203 IF ( ( zht_i_ini(i_hemis) > hi_max(jl-1) ) .AND. & 204 & ( zht_i_ini(i_hemis) <= hi_max(jl) ) ) THEN 205 jl0 = jl 134 135 !-------------------------------------------------------------------- 136 ! 3) Initialization of sea ice state variables 137 !-------------------------------------------------------------------- 138 IF( ln_limini_file )THEN 139 140 zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) 141 zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) 142 zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) 143 zts_u_ini(:,:) = si(jp_tsu)%fnow(:,:,1) 144 ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) 145 zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) 146 147 ELSE ! ln_limini_file = F 148 149 !----------------------------- 150 ! 3.1) Hemisphere-dependent arrays 151 !----------------------------- 152 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 153 DO jj = 1, jpj 154 DO ji = 1, jpi 155 IF( fcor(ji,jj) >= 0._wp ) THEN 156 zht_i_ini(ji,jj) = rn_hti_ini_n 157 zht_s_ini(ji,jj) = rn_hts_ini_n 158 zat_i_ini(ji,jj) = rn_ati_ini_n 159 zts_u_ini(ji,jj) = rn_tmi_ini_n 160 zsm_i_ini(ji,jj) = rn_smi_ini_n 161 ztm_i_ini(ji,jj) = rn_tmi_ini_n 162 ELSE 163 zht_i_ini(ji,jj) = rn_hti_ini_s 164 zht_s_ini(ji,jj) = rn_hts_ini_s 165 zat_i_ini(ji,jj) = rn_ati_ini_s 166 zts_u_ini(ji,jj) = rn_tmi_ini_s 167 zsm_i_ini(ji,jj) = rn_smi_ini_s 168 ztm_i_ini(ji,jj) = rn_tmi_ini_s 206 169 ENDIF 207 170 END DO 208 jl0 = MIN(jl0, i_fill) 209 210 !--- Concentrations 211 za_i_ini(jl0,i_hemis) = zat_i_ini(i_hemis) / SQRT(REAL(jpl)) 212 DO jl = 1, i_fill - 1 213 IF ( jl .NE. jl0 ) THEN 214 zsigma = 0.5 * zht_i_ini(i_hemis) 215 zarg = ( zh_i_ini(jl,i_hemis) - zht_i_ini(i_hemis) ) / zsigma 216 za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2) 217 ENDIF 218 END DO 219 220 zA = 0. ! sum of the areas in the jpl categories 221 DO jl = 1, i_fill - 1 222 zA = zA + za_i_ini(jl,i_hemis) 223 END DO 224 za_i_ini(i_fill,i_hemis) = zat_i_ini(i_hemis) - zA ! ice conc in the last category 225 IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 171 END DO 172 173 ENDIF ! ln_limini_file 226 174 227 !--- Ice thickness in the last category 228 zV = 0. ! sum of the volumes of the N-1 categories 229 DO jl = 1, i_fill - 1 230 zV = zV + za_i_ini(jl,i_hemis)*zh_i_ini(jl,i_hemis) 231 END DO 232 zh_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) 233 IF ( i_fill .LT. jpl ) zh_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 234 235 !--- volumes 236 zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zh_i_ini(:,i_hemis) 237 IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 238 239 ENDIF ! i_fill 240 241 !--------------------- 242 ! Compatibility tests 243 !--------------------- 244 ! Test 1: area conservation 245 zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons ) 246 IF ( zconv .LT. 1.0e-6 ) THEN 247 ztest_1 = 1 248 ELSE 249 ! this write is useful 250 IF(lwp) WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) 251 ztest_1 = 0 252 ENDIF 253 254 ! Test 2: volume conservation 255 zV_cons = SUM(zv_i_ini(:,i_hemis)) 256 zconv = ABS(zvt_i_ini(i_hemis) - zV_cons) 257 258 IF ( zconv .LT. 1.0e-6 ) THEN 259 ztest_2 = 1 260 ELSE 261 ! this write is useful 262 IF(lwp) WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 263 ' zvt_i_ini = ', zvt_i_ini(i_hemis) 264 ztest_2 = 0 265 ENDIF 266 267 ! Test 3: thickness of the last category is in-bounds ? 268 IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN 269 ztest_3 = 1 270 ELSE 271 ! this write is useful 272 IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(i_fill,i_hemis) = ', & 273 zh_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 274 ztest_3 = 0 275 ENDIF 276 277 ! Test 4: positivity of ice concentrations 278 ztest_4 = 1 279 DO jl = 1, jpl 280 IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 281 ! this write is useful 282 IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis) 283 ztest_4 = 0 284 ENDIF 285 END DO 286 287 ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 288 289 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 290 291 END DO ! i_fill 292 293 IF(lwp) THEN 294 WRITE(numout,*) ' ztests : ', ztests 295 IF ( ztests .NE. 4 ) THEN 296 WRITE(numout,*) 297 WRITE(numout,*) ' !!!! ALERT !!! ' 298 WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 299 WRITE(numout,*) 300 WRITE(numout,*) ' *** ztests is not equal to 4 ' 301 WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 302 WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(i_hemis) 303 WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(i_hemis) 304 ENDIF ! ztests .NE. 4 305 ENDIF 306 307 END DO ! i_hemis 308 309 !--------------------------------------------------------------------- 310 ! 3.3) Space-dependent arrays for ice state variables 311 !--------------------------------------------------------------------- 312 313 ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 314 DO jl = 1, jpl ! loop over categories 175 zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) ! ice volume 176 !--------------------------------------------------------------------- 177 ! 3.2) Distribute ice concentration and thickness into the categories 178 !--------------------------------------------------------------------- 179 ! a gaussian distribution for ice concentration is used 180 ! then we check whether the distribution fullfills 181 ! volume and area conservation, positivity and ice categories bounds 182 zh_i_ini(:,:,:) = 0._wp 183 za_i_ini(:,:,:) = 0._wp 184 zv_i_ini(:,:,:) = 0._wp 185 315 186 DO jj = 1, jpj 316 187 DO ji = 1, jpi 317 a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj)) ! concentration 318 ht_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj)) ! ice thickness 319 ht_s(ji,jj,jl) = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) ) ! snow depth 320 sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) ! salinity 321 o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp ! age (1 day) 322 t_su(ji,jj,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 323 324 ! This case below should not be used if (ht_s/ht_i) is ok in namelist 325 ! In case snow load is in excess that would lead to transformation from snow to ice 326 ! Then, transfer the snow excess into the ice (different from limthd_dh) 327 zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 328 ! recompute ht_i, ht_s avoiding out of bounds values 329 ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 330 ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 331 332 ! ice volume, salt content, age content 333 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume 334 v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume 335 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 336 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content 337 END DO 338 END DO 339 END DO 340 341 ! for constant salinity in time 342 IF( nn_icesal == 1 .OR. nn_icesal == 3 ) THEN 343 CALL lim_var_salprof 344 smv_i = sm_i * v_i 345 ENDIF 346 347 ! Snow temperature and heat content 348 DO jk = 1, nlay_s 188 189 IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 190 191 ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 192 ! ztests = 0 193 194 DO i_fill = jpl, 1, -1 195 196 ! IF( ztests .NE. 4 ) THEN 197 IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 198 !---------------------------- 199 ! fill the i_fill categories 200 !---------------------------- 201 ! *** 1 category to fill 202 IF ( i_fill .EQ. 1 ) THEN 203 zh_i_ini(ji,jj, 1) = zht_i_ini(ji,jj) 204 za_i_ini(ji,jj, 1) = zat_i_ini(ji,jj) 205 zh_i_ini(ji,jj,2:jpl) = 0._wp 206 za_i_ini(ji,jj,2:jpl) = 0._wp 207 ELSE 208 209 ! *** >1 categores to fill 210 !--- Ice thicknesses in the i_fill - 1 first categories 211 DO jl = 1, i_fill - 1 212 zh_i_ini(ji,jj,jl) = hi_mean(jl) 213 END DO 214 215 !--- jl0: most likely index where cc will be maximum 216 DO jl = 1, jpl 217 IF ( ( zht_i_ini(ji,jj) > hi_max(jl-1) ) .AND. & 218 & ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 219 jl0 = jl 220 ENDIF 221 END DO 222 jl0 = MIN(jl0, i_fill) 223 224 !--- Concentrations 225 za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 226 DO jl = 1, i_fill - 1 227 IF( jl .NE. jl0 )THEN 228 zsigma = 0.5 * zht_i_ini(ji,jj) 229 zarg = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / zsigma 230 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 231 ENDIF 232 END DO 233 234 zA = 0. ! sum of the areas in the jpl categories 235 DO jl = 1, i_fill - 1 236 zA = zA + za_i_ini(ji,jj,jl) 237 END DO 238 za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - zA ! ice conc in the last category 239 IF ( i_fill .LT. jpl ) za_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 240 241 !--- Ice thickness in the last category 242 zV = 0. ! sum of the volumes of the N-1 categories 243 DO jl = 1, i_fill - 1 244 zV = zV + za_i_ini(ji,jj,jl)*zh_i_ini(ji,jj,jl) 245 END DO 246 zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / za_i_ini(ji,jj,i_fill) 247 IF ( i_fill .LT. jpl ) zh_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 248 249 !--- volumes 250 zv_i_ini(ji,jj,:) = za_i_ini(ji,jj,:) * zh_i_ini(ji,jj,:) 251 IF ( i_fill .LT. jpl ) zv_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 252 253 ENDIF ! i_fill 254 255 !--------------------- 256 ! Compatibility tests 257 !--------------------- 258 ! Test 1: area conservation 259 zA_cons = SUM(za_i_ini(ji,jj,:)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons ) 260 IF ( zconv .LT. 1.0e-6 ) THEN 261 ztest_1 = 1 262 ELSE 263 !this write is useful 264 IF(lwp) WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(ji,jj) 265 ztest_1 = 0 266 ENDIF 267 268 ! Test 2: volume conservation 269 zV_cons = SUM(zv_i_ini(ji,jj,:)) 270 zconv = ABS(zvt_i_ini(ji,jj) - zV_cons) 271 272 IF( zconv .LT. 1.0e-6 ) THEN 273 ztest_2 = 1 274 ELSE 275 !this write is useful 276 IF(lwp) WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 277 ' zvt_i_ini = ', zvt_i_ini(ji,jj) 278 ztest_2 = 0 279 ENDIF 280 281 ! Test 3: thickness of the last category is in-bounds ? 282 IF ( zh_i_ini(ji,jj,i_fill) > hi_max(i_fill-1) ) THEN 283 ztest_3 = 1 284 ELSE 285 ! this write is useful 286 IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(ji,jj,i_fill) = ', & 287 zh_i_ini(ji,jj,i_fill), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 288 IF(lwp) WRITE(numout,*) ' ji,jj,i_fill ',ji,jj,i_fill 289 IF(lwp) WRITE(numout,*) 'zht_i_ini ',zht_i_ini(ji,jj) 290 ztest_3 = 0 291 ENDIF 292 293 ! Test 4: positivity of ice concentrations 294 ztest_4 = 1 295 DO jl = 1, jpl 296 IF ( za_i_ini(ji,jj,jl) .LT. 0._wp ) THEN 297 ! this write is useful 298 IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(ji,jj,jl) 299 ztest_4 = 0 300 ENDIF 301 END DO 302 303 ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 304 305 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 306 307 END DO ! i_fill 308 309 IF(lwp) THEN 310 WRITE(numout,*) ' ztests : ', ztests 311 IF( ztests .NE. 4 )THEN 312 WRITE(numout,*) 313 WRITE(numout,*) ' !!!! ALERT !!! ' 314 WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 315 WRITE(numout,*) 316 WRITE(numout,*) ' *** ztests is not equal to 4 ' 317 WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 318 WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 319 WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 320 ENDIF ! ztests .NE. 4 321 ENDIF 322 323 ENDIF ! zat_i_ini(ji,jj) > 0._wp .AND. zhm_i_ini(ji,jj) > 0._wp 324 325 ENDDO 326 ENDDO 327 328 !--------------------------------------------------------------------- 329 ! 3.3) Space-dependent arrays for ice state variables 330 !--------------------------------------------------------------------- 331 332 ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 349 333 DO jl = 1, jpl ! loop over categories 350 334 DO jj = 1, jpj 351 335 DO ji = 1, jpi 352 t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 353 ! Snow energy of melting 354 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 355 356 ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 357 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 336 a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini(ji,jj,jl) ! concentration 337 ht_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(ji,jj,jl) ! ice thickness 338 sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) ! salinity 339 o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp ! age (1 day) 340 t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 341 342 IF( zht_i_ini(ji,jj) > 0._wp )THEN 343 ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) ) ! snow depth 344 ELSE 345 ht_s(ji,jj,jl)= 0._wp 346 ENDIF 347 348 ! This case below should not be used if (ht_s/ht_i) is ok in namelist 349 ! In case snow load is in excess that would lead to transformation from snow to ice 350 ! Then, transfer the snow excess into the ice (different from limthd_dh) 351 zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 352 ! recompute ht_i, ht_s avoiding out of bounds values 353 ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 354 ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 355 356 ! ice volume, salt content, age content 357 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume 358 v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume 359 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 360 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content 358 361 END DO 359 362 END DO 360 363 END DO 361 END DO 362 363 ! Ice salinity, temperature and heat content 364 DO jk = 1, nlay_i 365 DO jl = 1, jpl ! loop over categories 366 DO jj = 1, jpj 367 DO ji = 1, jpi 368 t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 369 s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin 370 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 371 372 ! heat content per unit volume 373 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 374 + lfus * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 375 - rcp * ( ztmelts - rt0 ) ) 376 377 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 378 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 364 365 ! for constant salinity in time 366 IF( nn_icesal == 1 .OR. nn_icesal == 3 ) THEN 367 CALL lim_var_salprof 368 smv_i = sm_i * v_i 369 ENDIF 370 371 ! Snow temperature and heat content 372 DO jk = 1, nlay_s 373 DO jl = 1, jpl ! loop over categories 374 DO jj = 1, jpj 375 DO ji = 1, jpi 376 t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 377 ! Snow energy of melting 378 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 379 380 ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 381 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 382 END DO 379 383 END DO 380 384 END DO 381 385 END DO 382 END DO 383 384 tn_ice (:,:,:) = t_su (:,:,:) 385 386 ELSE 387 ! if ln_iceini=false 386 387 ! Ice salinity, temperature and heat content 388 DO jk = 1, nlay_i 389 DO jl = 1, jpl ! loop over categories 390 DO jj = 1, jpj 391 DO ji = 1, jpi 392 t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 393 s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin 394 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 395 396 ! heat content per unit volume 397 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 398 + lfus * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 399 - rcp * ( ztmelts - rt0 ) ) 400 401 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 402 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 403 END DO 404 END DO 405 END DO 406 END DO 407 408 tn_ice (:,:,:) = t_su (:,:,:) 409 410 ELSE ! if ln_limini=false 388 411 a_i (:,:,:) = 0._wp 389 412 v_i (:,:,:) = 0._wp … … 407 430 END DO 408 431 END DO 409 410 ENDIF ! ln_ iceini432 433 ENDIF ! ln_limini 411 434 412 435 at_i (:,:) = 0.0_wp … … 459 482 460 483 461 CALL wrk_dealloc( jpi, jpj, zswitch ) 462 CALL wrk_dealloc( jpi, jpj, zhemis ) 463 CALL wrk_dealloc( jpl, 2, zh_i_ini, za_i_ini, zv_i_ini ) 464 CALL wrk_dealloc( 2, zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 484 CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini, za_i_ini, zv_i_ini ) 485 CALL wrk_dealloc( jpi, jpj, zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 486 CALL wrk_dealloc( jpi, jpj, zswitch ) 465 487 466 488 END SUBROUTINE lim_istate … … 482 504 !! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization 483 505 !!----------------------------------------------------------------------------- 484 NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s, & 485 & rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s 486 INTEGER :: ios ! Local integer output status for namelist read 506 ! 507 INTEGER :: ios,ierr,inum_ice ! Local integer output status for namelist read 508 INTEGER :: ji,jj 509 INTEGER :: ifpr, ierror 510 ! 511 CHARACTER(len=100) :: cn_dir ! Root directory for location of ice files 512 TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi 513 TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read 514 ! 515 NAMELIST/namiceini/ ln_limini, ln_limini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, & 516 & rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, & 517 & rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s, & 518 & sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir 487 519 !!----------------------------------------------------------------------------- 488 520 ! … … 496 528 IF(lwm) WRITE ( numoni, namiceini ) 497 529 530 slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts 531 slf_i(jp_ati) = sn_ati ; slf_i(jp_tsu) = sn_tsu 532 slf_i(jp_tmi) = sn_tmi ; slf_i(jp_smi) = sn_smi 533 498 534 ! Define the initial parameters 499 535 ! ------------------------- … … 503 539 WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 504 540 WRITE(numout,*) '~~~~~~~~~~~~~~~' 505 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini 541 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_limini = ', ln_limini 542 WRITE(numout,*) ' ice initialization from a netcdf file ln_limini_file = ', ln_limini_file 506 543 WRITE(numout,*) ' threshold water temp. for initial sea-ice rn_thres_sst = ', rn_thres_sst 507 544 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n … … 517 554 ENDIF 518 555 556 IF( ln_limini_file ) THEN ! Ice initialization using input file 557 ! 558 ! set si structure 559 ALLOCATE( si(jpfldi), STAT=ierror ) 560 IF( ierror > 0 ) THEN 561 CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' ) ; RETURN 562 ENDIF 563 564 DO ifpr = 1, jpfldi 565 ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) ) 566 ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) 567 END DO 568 569 ! fill si with slf_i and control print 570 CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' ) 571 572 CALL fld_read( nit000, 1, si ) ! input fields provided at the current time-step 573 574 ENDIF 575 519 576 END SUBROUTINE lim_istate_init 520 577
Note: See TracChangeset
for help on using the changeset viewer.
