Changeset 5976
- Timestamp:
- 2015-12-02T12:46:21+01:00 (8 years ago)
- Location:
- branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM
- Files:
-
- 2 edited
-
CONFIG/SHARED/namelist_ice_lim3_ref (modified) (1 diff)
-
NEMO/LIM_SRC_3/limistate.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r5429 r5976 32 32 !------------------------------------------------------------------------------ 33 33 ln_iceini = .true. ! activate ice initialization (T) or not (F) 34 ln_iceini_file = .true. ! ice initialization from a netcdf file (T) or constant (F) 35 ! 36 ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 37 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 38 sn_hti = 'Ice_initialization' , -12 ,'hti' , .false. , .true. , 'yearly' , '' , '' , '' 39 sn_hts = 'Ice_initialization' , -12 ,'hts' , .false. , .true. , 'yearly' , '' , '' , '' 40 sn_ati = 'Ice_initialization' , -12 ,'ati' , .false. , .true. , 'yearly' , '' , '' , '' 41 sn_tsu = 'Ice_initialization' , -12 ,'tsu' , .false. , .true. , 'yearly' , '' , '' , '' 42 sn_tmi = 'Ice_initialization' , -12 ,'tmi' , .false. , .true. , 'yearly' , '' , '' , '' 43 sn_smi = 'Ice_initialization' , -12 ,'smi' , .false. , .true. , 'yearly' , '' , '' , '' 44 cn_dir='./' 45 ! 34 46 rn_thres_sst = 2.0 ! maximum water temperature with initial ice (degC) 35 47 rn_hts_ini_n = 0.3 ! initial real snow thickness (m), North -
branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r5541 r5976 28 28 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 29 29 USE wrk_nemo ! work arrays 30 USE fldread ! read input fields 31 USE iom 30 32 31 33 IMPLICIT NONE … … 47 49 REAL(wp) :: rn_tmi_ini_s ! initial temperature 48 50 49 LOGICAL :: ln_iceini ! initialization or not 51 INTEGER , PARAMETER :: jpfldi = 6 ! maximum number of files to read 52 INTEGER , PARAMETER :: jp_hti = 1 ! index of thick (m) at T-point 53 INTEGER , PARAMETER :: jp_hts = 2 ! index of thick (m) at T-point 54 INTEGER , PARAMETER :: jp_ati = 3 ! index of ice fraction (%) at T-point 55 INTEGER , PARAMETER :: jp_tsu = 4 ! index of ice surface temp (K) at T-point 56 INTEGER , PARAMETER :: jp_tmi = 5 ! index of ice temp at T-point 57 INTEGER , PARAMETER :: jp_smi = 6 ! index of ice sali at T-point 58 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read) 59 60 LOGICAL :: ln_iceini ! initialization or not 61 LOGICAL :: ln_iceini_file ! Ice initialization state from 2D netcdf file 50 62 !!---------------------------------------------------------------------- 51 63 !! LIM 3.0, UCL-LOCEAN-IPSL (2008) … … 88 100 REAL(wp) :: ztmelts, zdh 89 101 INTEGER :: i_hemis, i_fill, jl0 90 REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 91 REAL(wp), POINTER, DIMENSION(:) :: zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini 92 REAL(wp), POINTER, DIMENSION(:,:) :: zh_i_ini, za_i_ini, zv_i_ini 93 REAL(wp), POINTER, DIMENSION(:,:) :: zswitch ! ice indicator 94 INTEGER, POINTER, DIMENSION(:,:) :: zhemis ! hemispheric index 95 !-------------------------------------------------------------------- 96 97 CALL wrk_alloc( jpi, jpj, zswitch ) 98 CALL wrk_alloc( jpi, jpj, zhemis ) 99 CALL wrk_alloc( jpl, 2, zh_i_ini, za_i_ini, zv_i_ini ) 100 CALL wrk_alloc( 2, zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 102 REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 103 REAL(wp), POINTER, DIMENSION(:,:) :: zht_i_ini, zat_i_ini, zvt_i_ini !data from namelist or nc file 104 REAL(wp), POINTER, DIMENSION(:,:) :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 105 REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini, zv_i_ini !data by cattegories to fill 106 !-------------------------------------------------------------------- 107 108 CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini, za_i_ini, zv_i_ini ) 109 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 ) 101 110 102 111 IF(lwp) WRITE(numout,*) … … 117 126 118 127 ! basal temperature (considered at freezing point) 119 CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 120 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 128 t_bo(:,:) = ( eos_fzp( sss_m(:,:) ) + rt0 ) * tmask(:,:,1) 121 129 122 130 IF( ln_iceini ) THEN 123 131 124 !-------------------------------------------------------------------- 125 ! 2) Basal temperature, ice mask and hemispheric index 126 !-------------------------------------------------------------------- 127 128 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 129 DO ji = 1, jpi 130 IF( ( sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 131 zswitch(ji,jj) = 0._wp * tmask(ji,jj,1) ! no ice 132 ELSE 133 zswitch(ji,jj) = 1._wp * tmask(ji,jj,1) ! ice 134 ENDIF 132 !-------------------------------------------------------------------- 133 ! 2) Basal temperature, ice mask and hemispheric index 134 !-------------------------------------------------------------------- 135 136 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 137 DO ji = 1, jpi 138 IF( ( sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 139 zswitch(ji,jj) = 0._wp * tmask(ji,jj,1) ! no ice 140 ELSE 141 zswitch(ji,jj) = 1._wp * tmask(ji,jj,1) ! ice 142 ENDIF 143 END DO 135 144 END DO 136 END DO 137 138 139 ! Hemispheric index 140 DO jj = 1, jpj 141 DO ji = 1, jpi 142 IF( fcor(ji,jj) >= 0._wp ) THEN 143 zhemis(ji,jj) = 1 ! Northern hemisphere 144 ELSE 145 zhemis(ji,jj) = 2 ! Southern hemisphere 146 ENDIF 147 END DO 148 END DO 149 150 !-------------------------------------------------------------------- 151 ! 3) Initialization of sea ice state variables 152 !-------------------------------------------------------------------- 153 154 !----------------------------- 155 ! 3.1) Hemisphere-dependent arrays 156 !----------------------------- 157 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 158 zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s ! ice thickness 159 zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s ! snow depth 160 zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s ! ice concentration 161 zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s ! bulk ice salinity 162 ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s ! temperature (ice and snow) 163 164 zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:) ! ice volume 165 166 !--------------------------------------------------------------------- 167 ! 3.2) Distribute ice concentration and thickness into the categories 168 !--------------------------------------------------------------------- 169 ! a gaussian distribution for ice concentration is used 170 ! then we check whether the distribution fullfills 171 ! volume and area conservation, positivity and ice categories bounds 172 DO i_hemis = 1, 2 173 174 ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 175 176 ! note for the great nemo engineers: 177 ! only very few of the WRITE statements are necessary for the reference version 178 ! they were one day useful, but now i personally doubt of their 179 ! potential for bringing anything useful 180 181 DO i_fill = jpl, 1, -1 182 IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 183 !---------------------------- 184 ! fill the i_fill categories 185 !---------------------------- 186 ! *** 1 category to fill 187 IF ( i_fill .EQ. 1 ) THEN 188 zh_i_ini(1,i_hemis) = zht_i_ini(i_hemis) 189 za_i_ini(1,i_hemis) = zat_i_ini(i_hemis) 190 zh_i_ini(2:jpl,i_hemis) = 0._wp 191 za_i_ini(2:jpl,i_hemis) = 0._wp 192 ELSE 193 194 ! *** >1 categores to fill 195 !--- Ice thicknesses in the i_fill - 1 first categories 196 DO jl = 1, i_fill - 1 197 zh_i_ini(jl,i_hemis) = hi_mean(jl) 198 END DO 199 200 !--- jl0: most likely index where cc will be maximum 201 DO jl = 1, jpl 202 IF ( ( zht_i_ini(i_hemis) > hi_max(jl-1) ) .AND. & 203 & ( zht_i_ini(i_hemis) <= hi_max(jl) ) ) THEN 204 jl0 = jl 145 146 !-------------------------------------------------------------------- 147 ! 3) Initialization of sea ice state variables 148 !-------------------------------------------------------------------- 149 IF( ln_iceini_file )THEN 150 151 zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) 152 zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) 153 zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) 154 zts_u_ini(:,:) = si(jp_tsu)%fnow(:,:,1) 155 ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) 156 zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) 157 158 ELSE ! ln_iceini_file = F 159 160 !----------------------------- 161 ! 3.1) Hemisphere-dependent arrays 162 !----------------------------- 163 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 164 DO jj = 1, jpj 165 DO ji = 1, jpi 166 IF( fcor(ji,jj) >= 0._wp ) THEN 167 zht_i_ini(ji,jj) = rn_hti_ini_n 168 zht_s_ini(ji,jj) = rn_hts_ini_n 169 zat_i_ini(ji,jj) = rn_ati_ini_n 170 zsm_i_ini(ji,jj) = rn_smi_ini_n 171 ztm_i_ini(ji,jj) = rn_tmi_ini_n 172 ELSE 173 zht_i_ini(ji,jj) = rn_hti_ini_s 174 zht_s_ini(ji,jj) = rn_hts_ini_s 175 zat_i_ini(ji,jj) = rn_ati_ini_s 176 zsm_i_ini(ji,jj) = rn_smi_ini_s 177 ztm_i_ini(ji,jj) = rn_tmi_ini_s 205 178 ENDIF 206 179 END DO 207 jl0 = MIN(jl0, i_fill) 208 209 !--- Concentrations 210 za_i_ini(jl0,i_hemis) = zat_i_ini(i_hemis) / SQRT(REAL(jpl)) 211 DO jl = 1, i_fill - 1 212 IF ( jl .NE. jl0 ) THEN 213 zsigma = 0.5 * zht_i_ini(i_hemis) 214 zarg = ( zh_i_ini(jl,i_hemis) - zht_i_ini(i_hemis) ) / zsigma 215 za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2) 216 ENDIF 217 END DO 218 219 zA = 0. ! sum of the areas in the jpl categories 220 DO jl = 1, i_fill - 1 221 zA = zA + za_i_ini(jl,i_hemis) 222 END DO 223 za_i_ini(i_fill,i_hemis) = zat_i_ini(i_hemis) - zA ! ice conc in the last category 224 IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 225 226 !--- Ice thickness in the last category 227 zV = 0. ! sum of the volumes of the N-1 categories 228 DO jl = 1, i_fill - 1 229 zV = zV + za_i_ini(jl,i_hemis)*zh_i_ini(jl,i_hemis) 230 END DO 231 zh_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) 232 IF ( i_fill .LT. jpl ) zh_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 233 234 !--- volumes 235 zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zh_i_ini(:,i_hemis) 236 IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 237 238 ENDIF ! i_fill 239 240 !--------------------- 241 ! Compatibility tests 242 !--------------------- 243 ! Test 1: area conservation 244 zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons ) 245 IF ( zconv .LT. 1.0e-6 ) THEN 246 ztest_1 = 1 247 ELSE 248 ! this write is useful 249 IF(lwp) WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) 250 ztest_1 = 0 251 ENDIF 252 253 ! Test 2: volume conservation 254 zV_cons = SUM(zv_i_ini(:,i_hemis)) 255 zconv = ABS(zvt_i_ini(i_hemis) - zV_cons) 256 257 IF ( zconv .LT. 1.0e-6 ) THEN 258 ztest_2 = 1 259 ELSE 260 ! this write is useful 261 IF(lwp) WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 262 ' zvt_i_ini = ', zvt_i_ini(i_hemis) 263 ztest_2 = 0 264 ENDIF 265 266 ! Test 3: thickness of the last category is in-bounds ? 267 IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN 268 ztest_3 = 1 269 ELSE 270 ! this write is useful 271 IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(i_fill,i_hemis) = ', & 272 zh_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 273 ztest_3 = 0 274 ENDIF 275 276 ! Test 4: positivity of ice concentrations 277 ztest_4 = 1 278 DO jl = 1, jpl 279 IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 280 ! this write is useful 281 IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis) 282 ztest_4 = 0 283 ENDIF 284 END DO 285 286 ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 287 288 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 289 290 END DO ! i_fill 291 292 IF(lwp) THEN 293 WRITE(numout,*) ' ztests : ', ztests 294 IF ( ztests .NE. 4 ) THEN 295 WRITE(numout,*) 296 WRITE(numout,*) ' !!!! ALERT !!! ' 297 WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 298 WRITE(numout,*) 299 WRITE(numout,*) ' *** ztests is not equal to 4 ' 300 WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 301 WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(i_hemis) 302 WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(i_hemis) 303 ENDIF ! ztests .NE. 4 304 ENDIF 305 306 END DO ! i_hemis 307 308 !--------------------------------------------------------------------- 309 ! 3.3) Space-dependent arrays for ice state variables 310 !--------------------------------------------------------------------- 311 312 ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 313 DO jl = 1, jpl ! loop over categories 180 END DO 181 182 ENDIF ! ln_iceini_file 183 184 zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) ! ice volume 185 186 !--------------------------------------------------------------------- 187 ! 3.2) Distribute ice concentration and thickness into the categories 188 !--------------------------------------------------------------------- 189 ! a gaussian distribution for ice concentration is used 190 ! then we check whether the distribution fullfills 191 ! volume and area conservation, positivity and ice categories bounds 192 zh_i_ini(:,:,:) = 0._wp 193 za_i_ini(:,:,:) = 0._wp 194 zv_i_ini(:,:,:) = 0._wp 195 314 196 DO jj = 1, jpj 315 197 DO ji = 1, jpi 316 a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj)) ! concentration 317 ht_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj)) ! ice thickness 318 ht_s(ji,jj,jl) = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) ) ! snow depth 319 sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) ! salinity 320 o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp ! age (1 day) 321 t_su(ji,jj,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 322 323 ! This case below should not be used if (ht_s/ht_i) is ok in namelist 324 ! In case snow load is in excess that would lead to transformation from snow to ice 325 ! Then, transfer the snow excess into the ice (different from limthd_dh) 326 zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 327 ! recompute ht_i, ht_s avoiding out of bounds values 328 ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 329 ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 330 331 ! ice volume, salt content, age content 332 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume 333 v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume 334 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 335 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content 336 END DO 337 END DO 338 END DO 339 340 ! Snow temperature and heat content 341 DO jk = 1, nlay_s 198 199 IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 200 201 ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 202 ztests = 0 203 204 DO i_fill = jpl, 1, -1 205 206 IF( ztests .NE. 4 ) THEN 207 !---------------------------- 208 ! fill the i_fill categories 209 !---------------------------- 210 ! *** 1 category to fill 211 IF ( i_fill .EQ. 1 ) THEN 212 zh_i_ini(ji,jj, 1) = zht_i_ini(ji,jj) 213 za_i_ini(ji,jj, 1) = zat_i_ini(ji,jj) 214 zh_i_ini(ji,jj,2:jpl) = 0._wp 215 za_i_ini(ji,jj,2:jpl) = 0._wp 216 ELSE 217 218 ! *** >1 categores to fill 219 !--- Ice thicknesses in the i_fill - 1 first categories 220 DO jl = 1, i_fill - 1 221 zh_i_ini(ji,jj,jl) = hi_mean(jl) 222 END DO 223 224 !--- jl0: most likely index where cc will be maximum 225 DO jl = 1, jpl 226 IF ( ( zht_i_ini(ji,jj) > hi_max(jl-1) ) .AND. & 227 & ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 228 jl0 = jl 229 ENDIF 230 END DO 231 jl0 = MIN(jl0, i_fill) 232 233 !--- Concentrations 234 za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 235 DO jl = 1, i_fill - 1 236 IF( jl .NE. jl0 )THEN 237 zsigma = 0.5 * zht_i_ini(ji,jj) 238 zarg = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / zsigma 239 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 240 ENDIF 241 END DO 242 243 zA = 0. ! sum of the areas in the jpl categories 244 DO jl = 1, i_fill - 1 245 zA = zA + za_i_ini(ji,jj,jl) 246 END DO 247 za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - zA ! ice conc in the last category 248 IF ( i_fill .LT. jpl ) za_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 249 250 !--- Ice thickness in the last category 251 zV = 0. ! sum of the volumes of the N-1 categories 252 DO jl = 1, i_fill - 1 253 zV = zV + za_i_ini(ji,jj,jl)*zh_i_ini(ji,jj,jl) 254 END DO 255 zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / za_i_ini(ji,jj,i_fill) 256 IF ( i_fill .LT. jpl ) zh_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 257 258 !--- volumes 259 zv_i_ini(ji,jj,:) = za_i_ini(ji,jj,:) * zh_i_ini(ji,jj,:) 260 IF ( i_fill .LT. jpl ) zv_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 261 262 ENDIF ! i_fill 263 264 !--------------------- 265 ! Compatibility tests 266 !--------------------- 267 ! Test 1: area conservation 268 zA_cons = SUM(za_i_ini(ji,jj,:)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons ) 269 IF ( zconv .LT. 1.0e-6 ) THEN 270 ztest_1 = 1 271 ELSE 272 !this write is useful 273 IF(lwp) WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(ji,jj) 274 ztest_1 = 0 275 ENDIF 276 277 ! Test 2: volume conservation 278 zV_cons = SUM(zv_i_ini(ji,jj,:)) 279 zconv = ABS(zvt_i_ini(ji,jj) - zV_cons) 280 281 IF( zconv .LT. 1.0e-6 ) THEN 282 ztest_2 = 1 283 ELSE 284 !this write is useful 285 IF(lwp) WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 286 ' zvt_i_ini = ', zvt_i_ini(ji,jj) 287 ztest_2 = 0 288 ENDIF 289 290 ! Test 3: thickness of the last category is in-bounds ? 291 IF ( zh_i_ini(ji,jj,i_fill) > hi_max(i_fill-1) ) THEN 292 ztest_3 = 1 293 ELSE 294 ! this write is useful 295 IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(ji,jj,i_fill) = ', & 296 zh_i_ini(ji,jj,i_fill), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 297 ztest_3 = 0 298 ENDIF 299 300 ! Test 4: positivity of ice concentrations 301 ztest_4 = 1 302 DO jl = 1, jpl 303 IF ( za_i_ini(ji,jj,jl) .LT. 0._wp ) THEN 304 ! this write is useful 305 IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(ji,jj,jl) 306 ztest_4 = 0 307 ENDIF 308 END DO 309 310 ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 311 312 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 313 314 END DO ! i_fill 315 316 IF(lwp) THEN 317 WRITE(numout,*) ' ztests : ', ztests 318 IF( ztests .NE. 4 )THEN 319 WRITE(numout,*) 320 WRITE(numout,*) ' !!!! ALERT !!! ' 321 WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 322 WRITE(numout,*) 323 WRITE(numout,*) ' *** ztests is not equal to 4 ' 324 WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 325 WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 326 WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 327 ENDIF ! ztests .NE. 4 328 ENDIF 329 330 ENDIF ! zat_i_ini(ji,jj) > 0._wp .AND. zhm_i_ini(ji,jj) > 0._wp 331 332 ENDDO 333 ENDDO 334 335 !--------------------------------------------------------------------- 336 ! 3.3) Space-dependent arrays for ice state variables 337 !--------------------------------------------------------------------- 338 339 ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 342 340 DO jl = 1, jpl ! loop over categories 343 341 DO jj = 1, jpj 344 342 DO ji = 1, jpi 345 t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 346 ! Snow energy of melting 347 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 348 349 ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 350 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 343 a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini(ji,jj,jl) ! concentration 344 ht_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(ji,jj,jl) ! ice thickness 345 sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) ! salinity 346 o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp ! age (1 day) 347 t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 348 349 IF( zht_i_ini(ji,jj) > 0._wp )THEN 350 ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) ) ! snow depth 351 ELSE 352 ht_s(ji,jj,jl)= 0._wp 353 ENDIF 354 355 ! This case below should not be used if (ht_s/ht_i) is ok in namelist 356 ! In case snow load is in excess that would lead to transformation from snow to ice 357 ! Then, transfer the snow excess into the ice (different from limthd_dh) 358 zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 359 ! recompute ht_i, ht_s avoiding out of bounds values 360 ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 361 ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 362 363 ! ice volume, salt content, age content 364 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume 365 v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume 366 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 367 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content 351 368 END DO 352 369 END DO 353 370 END DO 354 END DO 355 356 ! Ice salinity, temperature and heat content 357 DO jk = 1, nlay_i 358 DO jl = 1, jpl ! loop over categories 359 DO jj = 1, jpj 360 DO ji = 1, jpi 361 t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 362 s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin 363 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 364 365 ! heat content per unit volume 366 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 367 + lfus * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 368 - rcp * ( ztmelts - rt0 ) ) 369 370 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 371 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 371 372 ! Snow temperature and heat content 373 DO jk = 1, nlay_s 374 DO jl = 1, jpl ! loop over categories 375 DO jj = 1, jpj 376 DO ji = 1, jpi 377 t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 378 ! Snow energy of melting 379 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 380 381 ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 382 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 383 END DO 372 384 END DO 373 385 END DO 374 386 END DO 375 END DO 376 377 tn_ice (:,:,:) = t_su (:,:,:) 378 379 ELSE 380 ! if ln_iceini=false 387 388 ! Ice salinity, temperature and heat content 389 DO jk = 1, nlay_i 390 DO jl = 1, jpl ! loop over categories 391 DO jj = 1, jpj 392 DO ji = 1, jpi 393 t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 394 s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin 395 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 396 397 ! heat content per unit volume 398 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 399 + lfus * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 400 - rcp * ( ztmelts - rt0 ) ) 401 402 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 403 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 404 END DO 405 END DO 406 END DO 407 END DO 408 409 tn_ice (:,:,:) = t_su (:,:,:) 410 411 ELSE ! if ln_iceini=false 381 412 a_i (:,:,:) = 0._wp 382 413 v_i (:,:,:) = 0._wp … … 400 431 END DO 401 432 END DO 402 433 403 434 ENDIF ! ln_iceini 404 435 … … 452 483 453 484 454 CALL wrk_dealloc( jpi, jpj, zswitch ) 455 CALL wrk_dealloc( jpi, jpj, zhemis ) 456 CALL wrk_dealloc( jpl, 2, zh_i_ini, za_i_ini, zv_i_ini ) 457 CALL wrk_dealloc( 2, zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 485 CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini, za_i_ini, zv_i_ini ) 486 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 ) 458 487 459 488 END SUBROUTINE lim_istate … … 475 504 !! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization 476 505 !!----------------------------------------------------------------------------- 477 NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s, & 478 & rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s 479 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_iceini, ln_iceini_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 480 519 !!----------------------------------------------------------------------------- 481 520 ! … … 489 528 IF(lwm) WRITE ( numoni, namiceini ) 490 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 491 534 ! Define the initial parameters 492 535 ! ------------------------- … … 497 540 WRITE(numout,*) '~~~~~~~~~~~~~~~' 498 541 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini 542 WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file 499 543 WRITE(numout,*) ' threshold water temp. for initial sea-ice rn_thres_sst = ', rn_thres_sst 500 544 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n … … 510 554 ENDIF 511 555 556 IF( ln_iceini_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 512 576 END SUBROUTINE lim_istate_init 513 577
Note: See TracChangeset
for help on using the changeset viewer.
