Changeset 6140 for trunk/NEMOGCM/NEMO/LIM_SRC_3
- Timestamp:
- 2015-12-21T12:35:23+01:00 (9 years ago)
- Location:
- trunk/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 1 deleted
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r5541 r6140 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 ice thickness (m) at T-point 53 INTEGER , PARAMETER :: jp_hts = 2 ! index of snow thicknes (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 102 REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 93 103 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 ) 104 REAL(wp), POINTER, DIMENSION(:,:) :: zht_i_ini, zat_i_ini, zvt_i_ini !data from namelist or nc file 105 REAL(wp), POINTER, DIMENSION(:,:) :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 106 REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini, zv_i_ini !data by cattegories to fill 107 !-------------------------------------------------------------------- 108 109 CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini, za_i_ini, zv_i_ini ) 110 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 ) 111 CALL wrk_alloc( jpi, jpj, zswitch ) 101 112 102 113 IF(lwp) WRITE(numout,*) … … 120 131 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 121 132 133 122 134 IF( ln_iceini ) THEN 123 135 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 136 !-------------------------------------------------------------------- 137 ! 2) Basal temperature, ice mask and hemispheric index 138 !-------------------------------------------------------------------- 139 140 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 141 DO ji = 1, jpi 142 IF( ( sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 143 zswitch(ji,jj) = 0._wp * tmask(ji,jj,1) ! no ice 144 ELSE 145 zswitch(ji,jj) = 1._wp * tmask(ji,jj,1) ! ice 146 ENDIF 147 END DO 135 148 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 149 150 !-------------------------------------------------------------------- 151 ! 3) Initialization of sea ice state variables 152 !-------------------------------------------------------------------- 153 IF( ln_iceini_file )THEN 154 155 zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) 156 zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) 157 zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) 158 zts_u_ini(:,:) = si(jp_tsu)%fnow(:,:,1) 159 ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) 160 zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) 161 162 ELSE ! ln_iceini_file = F 163 164 !----------------------------- 165 ! 3.1) Hemisphere-dependent arrays 166 !----------------------------- 167 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 168 DO jj = 1, jpj 169 DO ji = 1, jpi 170 IF( fcor(ji,jj) >= 0._wp ) THEN 171 zht_i_ini(ji,jj) = rn_hti_ini_n 172 zht_s_ini(ji,jj) = rn_hts_ini_n 173 zat_i_ini(ji,jj) = rn_ati_ini_n 174 zts_u_ini(ji,jj) = rn_tmi_ini_n 175 zsm_i_ini(ji,jj) = rn_smi_ini_n 176 ztm_i_ini(ji,jj) = rn_tmi_ini_n 177 ELSE 178 zht_i_ini(ji,jj) = rn_hti_ini_s 179 zht_s_ini(ji,jj) = rn_hts_ini_s 180 zat_i_ini(ji,jj) = rn_ati_ini_s 181 zts_u_ini(ji,jj) = rn_tmi_ini_s 182 zsm_i_ini(ji,jj) = rn_smi_ini_s 183 ztm_i_ini(ji,jj) = rn_tmi_ini_s 205 184 ENDIF 206 185 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 186 END DO 187 188 ENDIF ! ln_iceini_file 189 190 zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) ! ice volume 191 192 !--------------------------------------------------------------------- 193 ! 3.2) Distribute ice concentration and thickness into the categories 194 !--------------------------------------------------------------------- 195 ! a gaussian distribution for ice concentration is used 196 ! then we check whether the distribution fullfills 197 ! volume and area conservation, positivity and ice categories bounds 198 zh_i_ini(:,:,:) = 0._wp 199 za_i_ini(:,:,:) = 0._wp 200 zv_i_ini(:,:,:) = 0._wp 201 314 202 DO jj = 1, jpj 315 203 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 204 205 IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 206 207 ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 208 ! ztests = 0 209 210 DO i_fill = jpl, 1, -1 211 212 ! IF( ztests .NE. 4 ) THEN 213 IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 214 !---------------------------- 215 ! fill the i_fill categories 216 !---------------------------- 217 ! *** 1 category to fill 218 IF ( i_fill .EQ. 1 ) THEN 219 zh_i_ini(ji,jj, 1) = zht_i_ini(ji,jj) 220 za_i_ini(ji,jj, 1) = zat_i_ini(ji,jj) 221 zh_i_ini(ji,jj,2:jpl) = 0._wp 222 za_i_ini(ji,jj,2:jpl) = 0._wp 223 ELSE 224 225 ! *** >1 categores to fill 226 !--- Ice thicknesses in the i_fill - 1 first categories 227 DO jl = 1, i_fill - 1 228 zh_i_ini(ji,jj,jl) = hi_mean(jl) 229 END DO 230 231 !--- jl0: most likely index where cc will be maximum 232 DO jl = 1, jpl 233 IF ( ( zht_i_ini(ji,jj) > hi_max(jl-1) ) .AND. & 234 & ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 235 jl0 = jl 236 ENDIF 237 END DO 238 jl0 = MIN(jl0, i_fill) 239 240 !--- Concentrations 241 za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 242 DO jl = 1, i_fill - 1 243 IF( jl .NE. jl0 )THEN 244 zsigma = 0.5 * zht_i_ini(ji,jj) 245 zarg = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / zsigma 246 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 247 ENDIF 248 END DO 249 250 zA = 0. ! sum of the areas in the jpl categories 251 DO jl = 1, i_fill - 1 252 zA = zA + za_i_ini(ji,jj,jl) 253 END DO 254 za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - zA ! ice conc in the last category 255 IF ( i_fill .LT. jpl ) za_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 256 257 !--- Ice thickness in the last category 258 zV = 0. ! sum of the volumes of the N-1 categories 259 DO jl = 1, i_fill - 1 260 zV = zV + za_i_ini(ji,jj,jl)*zh_i_ini(ji,jj,jl) 261 END DO 262 zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / za_i_ini(ji,jj,i_fill) 263 IF ( i_fill .LT. jpl ) zh_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 264 265 !--- volumes 266 zv_i_ini(ji,jj,:) = za_i_ini(ji,jj,:) * zh_i_ini(ji,jj,:) 267 IF ( i_fill .LT. jpl ) zv_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 268 269 ENDIF ! i_fill 270 271 !--------------------- 272 ! Compatibility tests 273 !--------------------- 274 ! Test 1: area conservation 275 zA_cons = SUM(za_i_ini(ji,jj,:)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons ) 276 IF ( zconv .LT. 1.0e-6 ) THEN 277 ztest_1 = 1 278 ELSE 279 !this write is useful 280 IF(lwp) WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(ji,jj) 281 ztest_1 = 0 282 ENDIF 283 284 ! Test 2: volume conservation 285 zV_cons = SUM(zv_i_ini(ji,jj,:)) 286 zconv = ABS(zvt_i_ini(ji,jj) - zV_cons) 287 288 IF( zconv .LT. 1.0e-6 ) THEN 289 ztest_2 = 1 290 ELSE 291 !this write is useful 292 IF(lwp) WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 293 ' zvt_i_ini = ', zvt_i_ini(ji,jj) 294 ztest_2 = 0 295 ENDIF 296 297 ! Test 3: thickness of the last category is in-bounds ? 298 IF ( zh_i_ini(ji,jj,i_fill) > hi_max(i_fill-1) ) THEN 299 ztest_3 = 1 300 ELSE 301 ! this write is useful 302 IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(ji,jj,i_fill) = ', & 303 zh_i_ini(ji,jj,i_fill), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 304 IF(lwp) WRITE(numout,*) ' ji,jj,i_fill ',ji,jj,i_fill 305 IF(lwp) WRITE(numout,*) 'zht_i_ini ',zht_i_ini(ji,jj) 306 ztest_3 = 0 307 ENDIF 308 309 ! Test 4: positivity of ice concentrations 310 ztest_4 = 1 311 DO jl = 1, jpl 312 IF ( za_i_ini(ji,jj,jl) .LT. 0._wp ) THEN 313 ! this write is useful 314 IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(ji,jj,jl) 315 ztest_4 = 0 316 ENDIF 317 END DO 318 319 ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 320 321 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 322 323 END DO ! i_fill 324 325 IF(lwp) THEN 326 WRITE(numout,*) ' ztests : ', ztests 327 IF( ztests .NE. 4 )THEN 328 WRITE(numout,*) 329 WRITE(numout,*) ' !!!! ALERT !!! ' 330 WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 331 WRITE(numout,*) 332 WRITE(numout,*) ' *** ztests is not equal to 4 ' 333 WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 334 WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 335 WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 336 ENDIF ! ztests .NE. 4 337 ENDIF 338 339 ENDIF ! zat_i_ini(ji,jj) > 0._wp .AND. zhm_i_ini(ji,jj) > 0._wp 340 341 ENDDO 342 ENDDO 343 344 !--------------------------------------------------------------------- 345 ! 3.3) Space-dependent arrays for ice state variables 346 !--------------------------------------------------------------------- 347 348 ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 342 349 DO jl = 1, jpl ! loop over categories 343 350 DO jj = 1, jpj 344 351 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 352 a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini(ji,jj,jl) ! concentration 353 ht_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(ji,jj,jl) ! ice thickness 354 sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) ! salinity 355 o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp ! age (1 day) 356 t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 357 358 IF( zht_i_ini(ji,jj) > 0._wp )THEN 359 ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) ) ! snow depth 360 ELSE 361 ht_s(ji,jj,jl)= 0._wp 362 ENDIF 363 364 ! This case below should not be used if (ht_s/ht_i) is ok in namelist 365 ! In case snow load is in excess that would lead to transformation from snow to ice 366 ! Then, transfer the snow excess into the ice (different from limthd_dh) 367 zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 368 ! recompute ht_i, ht_s avoiding out of bounds values 369 ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 370 ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 371 372 ! ice volume, salt content, age content 373 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume 374 v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume 375 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 376 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content 351 377 END DO 352 378 END DO 353 379 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 380 381 ! Snow temperature and heat content 382 DO jk = 1, nlay_s 383 DO jl = 1, jpl ! loop over categories 384 DO jj = 1, jpj 385 DO ji = 1, jpi 386 t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 387 ! Snow energy of melting 388 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 389 390 ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 391 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 392 END DO 372 393 END DO 373 394 END DO 374 395 END DO 375 END DO 376 377 tn_ice (:,:,:) = t_su (:,:,:) 378 379 ELSE 380 ! if ln_iceini=false 396 397 ! Ice salinity, temperature and heat content 398 DO jk = 1, nlay_i 399 DO jl = 1, jpl ! loop over categories 400 DO jj = 1, jpj 401 DO ji = 1, jpi 402 t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 403 s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin 404 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 405 406 ! heat content per unit volume 407 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 408 + lfus * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 409 - rcp * ( ztmelts - rt0 ) ) 410 411 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 412 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 413 END DO 414 END DO 415 END DO 416 END DO 417 418 tn_ice (:,:,:) = t_su (:,:,:) 419 420 ELSE ! if ln_iceini=false 381 421 a_i (:,:,:) = 0._wp 382 422 v_i (:,:,:) = 0._wp … … 400 440 END DO 401 441 END DO 402 442 403 443 ENDIF ! ln_iceini 404 444 … … 452 492 453 493 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 ) 494 CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini, za_i_ini, zv_i_ini ) 495 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 ) 496 CALL wrk_dealloc( jpi, jpj, zswitch ) 458 497 459 498 END SUBROUTINE lim_istate … … 475 514 !! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization 476 515 !!----------------------------------------------------------------------------- 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 516 ! 517 INTEGER :: ios,ierr,inum_ice ! Local integer output status for namelist read 518 INTEGER :: ji,jj 519 INTEGER :: ifpr, ierror 520 ! 521 CHARACTER(len=100) :: cn_dir ! Root directory for location of ice files 522 TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi 523 TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read 524 ! 525 NAMELIST/namiceini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, & 526 & rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, & 527 & rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s, & 528 & sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir 480 529 !!----------------------------------------------------------------------------- 481 530 ! … … 489 538 IF(lwm) WRITE ( numoni, namiceini ) 490 539 540 slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts 541 slf_i(jp_ati) = sn_ati ; slf_i(jp_tsu) = sn_tsu 542 slf_i(jp_tmi) = sn_tmi ; slf_i(jp_smi) = sn_smi 543 491 544 ! Define the initial parameters 492 545 ! ------------------------- … … 497 550 WRITE(numout,*) '~~~~~~~~~~~~~~~' 498 551 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini 552 WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file 499 553 WRITE(numout,*) ' threshold water temp. for initial sea-ice rn_thres_sst = ', rn_thres_sst 500 554 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n … … 510 564 ENDIF 511 565 566 IF( ln_iceini_file ) THEN ! Ice initialization using input file 567 ! 568 ! set si structure 569 ALLOCATE( si(jpfldi), STAT=ierror ) 570 IF( ierror > 0 ) THEN 571 CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' ) ; RETURN 572 ENDIF 573 574 DO ifpr = 1, jpfldi 575 ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) ) 576 ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) 577 END DO 578 579 ! fill si with slf_i and control print 580 CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' ) 581 582 CALL fld_read( nit000, 1, si ) ! input fields provided at the current time-step 583 584 ENDIF 585 512 586 END SUBROUTINE lim_istate_init 513 587 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
r5341 r6140 77 77 WRITE(numout,*) 78 78 SELECT CASE ( jprstlib ) 79 CASE ( jprstdimg )80 WRITE(numout,*) ' open ice restart binary file: ',TRIM(clpath)//clname81 79 CASE DEFAULT 82 80 WRITE(numout,*) ' open ice restart NetCDF file: ',TRIM(clpath)//clname … … 331 329 ENDIF 332 330 333 IF ( jprstlib == jprstdimg ) THEN334 ! eventually read netcdf file (monobloc) for restarting on different number of processors335 ! if {cn_icerst_in}.nc exists, then set jlibalt to jpnf90336 INQUIRE( FILE = TRIM(cn_icerst_indir)//'/'//TRIM(cn_icerst_in)//'.nc', EXIST = llok )337 IF ( llok ) THEN ; jlibalt = jpnf90 ; ELSE ; jlibalt = jprstlib ; ENDIF338 ENDIF339 340 331 CALL iom_open ( TRIM(cn_icerst_indir)//'/'//cn_icerst_in, numrir, kiolib = jprstlib ) 341 332 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r5407 r6140 23 23 !! lim_sbc_tau : update i- and j-stresses, and its modulus at the ocean surface 24 24 !!---------------------------------------------------------------------- 25 USE par_oce ! ocean parameters 26 USE phycst ! physical constants 27 USE dom_oce ! ocean domain 28 USE ice ! LIM sea-ice variables 29 USE sbc_ice ! Surface boundary condition: sea-ice fields 30 USE sbc_oce ! Surface boundary condition: ocean fields 31 USE sbccpl 32 USE oce , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 33 USE albedo ! albedo parameters 34 USE lbclnk ! ocean lateral boundary condition - MPP exchanges 35 USE lib_mpp ! MPP library 36 USE wrk_nemo ! work arrays 37 USE in_out_manager ! I/O manager 38 USE prtctl ! Print control 39 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 40 USE traqsr ! add penetration of solar flux in the calculation of heat budget 41 USE iom 42 USE domvvl ! Variable volume 43 USE limctl 44 USE limcons 25 USE par_oce ! ocean parameters 26 USE oce , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 27 USE phycst ! physical constants 28 USE dom_oce ! ocean domain 29 USE ice ! LIM sea-ice variables 30 USE sbc_ice ! Surface boundary condition: sea-ice fields 31 USE sbc_oce ! Surface boundary condition: ocean fields 32 USE sbccpl ! Surface boundary condition: coupled interface 33 USE albedo ! albedo parameters 34 USE traqsr ! add penetration of solar flux in the calculation of heat budget 35 USE domvvl ! Variable volume 36 USE limctl ! 37 USE limcons ! 38 ! 39 USE in_out_manager ! I/O manager 40 USE iom ! xIO server 41 USE lbclnk ! ocean lateral boundary condition - MPP exchanges 42 USE lib_mpp ! MPP library 43 USE wrk_nemo ! work arrays 44 USE prtctl ! Print control 45 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 45 46 46 47 IMPLICIT NONE 47 48 PRIVATE 48 49 49 PUBLIC lim_sbc_init ! called by sbc _lim_init50 PUBLIC lim_sbc_init ! called by sbcice_lim 50 51 PUBLIC lim_sbc_flx ! called by sbc_ice_lim 51 52 PUBLIC lim_sbc_tau ! called by sbc_ice_lim … … 57 58 !! * Substitutions 58 59 # include "vectopt_loop_substitute.h90" 59 # include "domzgr_substitute.h90"60 60 !!---------------------------------------------------------------------- 61 61 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 101 101 !! The ref should be Rousset et al., 2015 102 102 !!--------------------------------------------------------------------- 103 INTEGER, INTENT(in) :: kt 104 INTEGER :: ji, jj, jl, jk ! dummy loop indices105 REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2)106 REAL(wp) :: zq sr ! New solar flux received by the ocean107 !108 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_cs, zalb_os ! 2D/3D workspace103 INTEGER, INTENT(in) :: kt ! number of iteration 104 ! 105 INTEGER :: ji, jj, jl, jk ! dummy loop indices 106 REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2) 107 REAL(wp) :: zqsr ! New solar flux received by the ocean 108 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_cs, zalb_os ! 3D workspace 109 109 !!--------------------------------------------------------------------- 110 110 ! 111 111 ! make calls for heat fluxes before it is modified 112 112 IF( iom_use('qsr_oce') ) CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) ) ! solar flux at ocean surface … … 198 198 ! Snow/ice albedo (only if sent to coupler, useless in forced mode) ! 199 199 !------------------------------------------------------------------------! 200 CALL wrk_alloc( jpi, jpj, jpl,zalb_cs, zalb_os )200 CALL wrk_alloc( jpi,jpj,jpl, zalb_cs, zalb_os ) 201 201 CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 202 202 alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 203 CALL wrk_dealloc( jpi, jpj, jpl,zalb_cs, zalb_os )203 CALL wrk_dealloc( jpi,jpj,jpl, zalb_cs, zalb_os ) 204 204 205 205 ! conservation test 206 IF( ln_limdiahsb ) CALL lim_cons_final( 'limsbc' )206 IF( ln_limdiahsb ) CALL lim_cons_final( 'limsbc' ) 207 207 208 208 ! control prints 209 209 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 3, ' - Final state lim_sbc - ' ) 210 210 ! 211 211 IF(ln_ctl) THEN 212 212 CALL prt_ctl( tab2d_1=qsr , clinfo1=' lim_sbc: qsr : ', tab2d_2=qns , clinfo2=' qns : ' ) … … 215 215 CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl ) 216 216 ENDIF 217 217 ! 218 218 END SUBROUTINE lim_sbc_flx 219 219 … … 246 246 INTEGER , INTENT(in) :: kt ! ocean time-step index 247 247 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents 248 ! !248 ! 249 249 INTEGER :: ji, jj ! dummy loop indices 250 250 REAL(wp) :: zat_u, zutau_ice, zu_t, zmodt ! local scalar … … 303 303 !! ** input : Namelist namicedia 304 304 !!------------------------------------------------------------------- 305 INTEGER :: ji, jj, jk ! dummy loop indices 306 REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar 305 INTEGER :: ji, jj, jk ! dummy loop indices 306 REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar 307 !!------------------------------------------------------------------- 308 ! 307 309 IF(lwp) WRITE(numout,*) 308 310 IF(lwp) WRITE(numout,*) 'lim_sbc_init : LIM-3 sea-ice - surface boundary condition' … … 335 337 sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 336 338 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 337 #if defined key_vvl 338 ! key_vvl necessary? clem: yes for compilation purpose 339 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 340 fse3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 341 fse3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 342 ENDDO 343 fse3t_a(:,:,:) = fse3t_b(:,:,:) 344 ! Reconstruction of all vertical scale factors at now and before time 345 ! steps 346 ! ============================================================================= 347 ! Horizontal scale factor interpolations 348 ! -------------------------------------- 349 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 350 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 351 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 352 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 353 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 354 ! Vertical scale factor interpolations 355 ! ------------------------------------ 356 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 357 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 358 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 359 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 360 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 361 ! t- and w- points depth 362 ! ---------------------- 363 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 364 fsdepw_n(:,:,1) = 0.0_wp 365 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 366 DO jk = 2, jpk 367 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 368 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 369 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 370 END DO 371 #endif 339 340 !!gm I really don't like this staff here... Find a way to put that elsewhere or differently 341 !!gm 342 IF( .NOT.ln_linssh ) THEN 343 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 344 e3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 345 e3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 346 END DO 347 e3t_a(:,:,:) = e3t_b(:,:,:) 348 ! Reconstruction of all vertical scale factors at now and before time-steps 349 ! ========================================================================= 350 ! Horizontal scale factor interpolations 351 ! -------------------------------------- 352 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 353 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 354 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 355 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 356 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 357 ! Vertical scale factor interpolations 358 ! ------------------------------------ 359 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 360 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 361 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 362 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 363 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 364 ! t- and w- points depth 365 ! ---------------------- 366 !!gm not sure of that.... 367 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 368 gdepw_n(:,:,1) = 0.0_wp 369 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 370 DO jk = 2, jpk 371 gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk) 372 gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1) 373 gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:) 374 END DO 375 ENDIF 372 376 ENDIF 373 377 ENDIF ! .NOT. ln_rstart 374 378 ! 375 376 379 END SUBROUTINE lim_sbc_init 377 380 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r5836 r6140 25 25 USE sbc_oce ! Surface boundary condition: ocean fields 26 26 USE sbc_ice ! Surface boundary condition: ice fields 27 USE thd_ice ! LIM thermodynamic sea-ice variables28 USE dom_ice ! LIM sea-ice domain27 USE dom_ice ! LIM: sea-ice domain 28 USE thd_ice ! LIM: thermodynamic sea-ice variables 29 29 USE limthd_dif ! LIM: thermodynamics, vertical diffusion 30 30 USE limthd_dh ! LIM: thermodynamics, ice and snow thickness variation 31 31 USE limthd_sal ! LIM: thermodynamics, ice salinity 32 32 USE limthd_ent ! LIM: thermodynamics, ice enthalpy redistribution 33 USE limthd_lac ! LIM -3lateral accretion34 USE limitd_th ! remapping thickness distribution33 USE limthd_lac ! LIM: lateral accretion 34 USE limitd_th ! LIM: remapping thickness distribution 35 35 USE limtab ! LIM: 1D <==> 2D transformation 36 36 USE limvar ! LIM: sea-ice variables 37 USE limcons ! LIM: conservation tests 38 USE limctl ! LIM: control print 39 ! 40 USE in_out_manager ! I/O manager 41 USE prtctl ! Print control 37 42 USE lbclnk ! lateral boundary condition - MPP links 38 43 USE lib_mpp ! MPP library 39 44 USE wrk_nemo ! work arrays 40 USE in_out_manager ! I/O manager41 USE prtctl ! Print control42 45 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 43 46 USE timing ! Timing 44 USE limcons ! conservation tests45 USE limctl46 47 47 48 IMPLICIT NONE … … 52 53 53 54 !! * Substitutions 54 # include "domzgr_substitute.h90"55 55 # include "vectopt_loop_substitute.h90" 56 56 !!---------------------------------------------------------------------- … … 81 81 !!--------------------------------------------------------------------- 82 82 INTEGER, INTENT(in) :: kt ! number of iteration 83 ! !83 ! 84 84 INTEGER :: ji, jj, jk, jl ! dummy loop indices 85 85 INTEGER :: nbpb ! nb of icy pts for vertical thermo calculations 86 INTEGER :: ii, ij ! temporary dummy loop index87 86 REAL(wp) :: zfric_u, zqld, zqfr 88 87 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 89 88 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 90 89 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 91 !92 90 !!------------------------------------------------------------------- 93 91 94 IF( nn_timing == 1 ) CALL timing_start('limthd')92 IF( nn_timing == 1 ) CALL timing_start('limthd') 95 93 96 94 ! conservation test 97 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)95 IF( ln_limdiahsb ) CALL lim_cons_hsm( 0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 98 96 99 97 CALL lim_var_glo2eqv … … 147 145 148 146 ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 149 zqfr = tmask(ji,jj,1) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )147 zqfr = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 150 148 151 149 ! --- Energy from the turbulent oceanic heat flux (W/m2) --- ! … … 226 224 227 225 IF( nbpb > 0 ) THEN ! If there is no ice, do nothing. 228 229 !-------------------------! 230 ! --- Move to 1D arrays --- 231 !-------------------------! 232 CALL lim_thd_1d2d( nbpb, jl, 1 ) 233 234 !--------------------------------------! 235 ! --- Ice/Snow Temperature profile --- ! 236 !--------------------------------------! 237 CALL lim_thd_dif( 1, nbpb ) 238 239 !---------------------------------! 240 ! --- Ice/Snow thickness --- ! 241 !---------------------------------! 242 CALL lim_thd_dh( 1, nbpb ) 243 244 ! --- Ice enthalpy remapping --- ! 245 CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) 246 247 !---------------------------------! 248 ! --- Ice salinity --- ! 249 !---------------------------------! 250 CALL lim_thd_sal( 1, nbpb ) 251 252 !---------------------------------! 253 ! --- temperature update --- ! 254 !---------------------------------! 255 CALL lim_thd_temp( 1, nbpb ) 256 257 !------------------------------------! 258 ! --- lateral melting if monocat --- ! 259 !------------------------------------! 226 ! 227 CALL lim_thd_1d2d( nbpb, jl, 1 ) ! --- Move to 1D arrays ---! 228 ! 229 CALL lim_thd_dif ( 1, nbpb ) ! --- Ice/Snow Temperature profile --- ! 230 ! 231 CALL lim_thd_dh ( 1, nbpb ) ! --- Ice/Snow thickness ---! 232 ! 233 CALL lim_thd_ent ( 1, nbpb, q_i_1d(1:nbpb,:) ) ! --- Ice enthalpy remapping --- ! 234 ! 235 CALL lim_thd_sal ( 1, nbpb ) ! --- Ice salinity --- ! 236 ! 237 CALL lim_thd_temp( 1, nbpb ) ! --- temperature update --- ! 238 ! 239 ! ! --- lateral melting if monocat --- ! 260 240 IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN 261 241 CALL lim_thd_lam( 1, nbpb ) 262 242 END IF 263 264 !-------------------------! 265 ! --- Move to 2D arrays --- 266 !-------------------------! 267 CALL lim_thd_1d2d( nbpb, jl, 2 ) 268 269 ! 270 IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? 243 ! 244 CALL lim_thd_1d2d( nbpb, jl, 2 ) ! --- Move to 2D arrays --- 245 ! 246 IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? 271 247 ENDIF 272 248 ! … … 410 386 ENDIF 411 387 ! 412 IF( nn_timing == 1 ) CALL timing_stop('limthd')413 388 IF( nn_timing == 1 ) CALL timing_stop('limthd') 389 ! 414 390 END SUBROUTINE lim_thd 415 391 … … 424 400 !!------------------------------------------------------------------- 425 401 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 426 ! !402 ! 427 403 INTEGER :: ji, jk ! dummy loop indices 428 404 REAL(wp) :: ztmelts, zaaa, zbbb, zccc, zdiscrim ! local scalar … … 444 420 END DO 445 421 END DO 446 422 ! 447 423 END SUBROUTINE lim_thd_temp 424 448 425 449 426 SUBROUTINE lim_thd_lam( kideb, kiut ) … … 455 432 !!----------------------------------------------------------------------- 456 433 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 457 INTEGER :: ji ! dummy loop indices 458 REAL(wp) :: zhi_bef ! ice thickness before thermo 459 REAL(wp) :: zdh_mel, zda_mel ! net melting 460 REAL(wp) :: zvi, zvs ! ice/snow volumes 461 434 ! 435 INTEGER :: ji ! dummy loop indices 436 REAL(wp) :: zhi_bef ! ice thickness before thermo 437 REAL(wp) :: zdh_mel, zda_mel ! net melting 438 REAL(wp) :: zvi, zvs ! ice/snow volumes 439 !!----------------------------------------------------------------------- 440 ! 462 441 DO ji = kideb, kiut 463 442 zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) … … 477 456 END IF 478 457 END DO 479 458 ! 480 459 END SUBROUTINE lim_thd_lam 460 481 461 482 462 SUBROUTINE lim_thd_1d2d( nbpb, jl, kn ) … … 486 466 !! ** Purpose : move arrays from 1d to 2d and the reverse 487 467 !!----------------------------------------------------------------------- 488 INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D 489 ! 2= from 1D to 2D 468 INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D 490 469 INTEGER, INTENT(in) :: nbpb ! size of 1D arrays 491 470 INTEGER, INTENT(in) :: jl ! ice cat 471 ! 492 472 INTEGER :: jk ! dummy loop indices 493 473 !!----------------------------------------------------------------------- 474 ! 494 475 SELECT CASE( kn ) 495 496 CASE( 1 ) 497 476 ! 477 CASE( 1 ) ! from 2D to 1D 478 ! 498 479 CALL tab_2d_1d( nbpb, at_i_1d (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) ) 499 480 CALL tab_2d_1d( nbpb, a_i_1d (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 500 481 CALL tab_2d_1d( nbpb, ht_i_1d (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 501 482 CALL tab_2d_1d( nbpb, ht_s_1d (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 502 483 ! 503 484 CALL tab_2d_1d( nbpb, t_su_1d (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 504 485 CALL tab_2d_1d( nbpb, sm_i_1d (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) … … 512 493 CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 513 494 END DO 514 495 ! 515 496 CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) ) 516 497 CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) … … 526 507 CALL tab_2d_1d( nbpb, qlead_1d (1:nbpb), qlead , jpi, jpj, npb(1:nbpb) ) 527 508 CALL tab_2d_1d( nbpb, fhld_1d (1:nbpb), fhld , jpi, jpj, npb(1:nbpb) ) 528 509 ! 529 510 CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw , jpi, jpj, npb(1:nbpb) ) 530 511 CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub , jpi, jpj, npb(1:nbpb) ) 531 512 ! 532 513 CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog , jpi, jpj, npb(1:nbpb) ) 533 514 CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom , jpi, jpj, npb(1:nbpb) ) … … 536 517 CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res , jpi, jpj, npb(1:nbpb) ) 537 518 CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr , jpi, jpj, npb(1:nbpb) ) 538 519 ! 539 520 CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog , jpi, jpj, npb(1:nbpb) ) 540 521 CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom , jpi, jpj, npb(1:nbpb) ) … … 543 524 CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) ) 544 525 CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res , jpi, jpj, npb(1:nbpb) ) 545 526 ! 546 527 CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd , jpi, jpj, npb(1:nbpb) ) 547 528 CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr , jpi, jpj, npb(1:nbpb) ) … … 557 538 CALL tab_2d_1d( nbpb, hfx_err_dif_1d (1:nbpb), hfx_err_dif , jpi, jpj, npb(1:nbpb) ) 558 539 CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) ) 559 560 CASE( 2 ) 561 540 ! 541 CASE( 2 ) ! from 1D to 2D 542 ! 562 543 CALL tab_1d_2d( nbpb, at_i , npb, at_i_1d (1:nbpb) , jpi, jpj ) 563 544 CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_1d (1:nbpb) , jpi, jpj ) … … 576 557 END DO 577 558 CALL tab_1d_2d( nbpb, qlead , npb, qlead_1d (1:nbpb) , jpi, jpj ) 578 559 ! 579 560 CALL tab_1d_2d( nbpb, wfx_snw , npb, wfx_snw_1d(1:nbpb) , jpi, jpj ) 580 561 CALL tab_1d_2d( nbpb, wfx_sub , npb, wfx_sub_1d(1:nbpb) , jpi, jpj ) 581 562 ! 582 563 CALL tab_1d_2d( nbpb, wfx_bog , npb, wfx_bog_1d(1:nbpb) , jpi, jpj ) 583 564 CALL tab_1d_2d( nbpb, wfx_bom , npb, wfx_bom_1d(1:nbpb) , jpi, jpj ) … … 586 567 CALL tab_1d_2d( nbpb, wfx_res , npb, wfx_res_1d(1:nbpb) , jpi, jpj ) 587 568 CALL tab_1d_2d( nbpb, wfx_spr , npb, wfx_spr_1d(1:nbpb) , jpi, jpj ) 588 569 ! 589 570 CALL tab_1d_2d( nbpb, sfx_bog , npb, sfx_bog_1d(1:nbpb) , jpi, jpj ) 590 571 CALL tab_1d_2d( nbpb, sfx_bom , npb, sfx_bom_1d(1:nbpb) , jpi, jpj ) … … 593 574 CALL tab_1d_2d( nbpb, sfx_res , npb, sfx_res_1d(1:nbpb) , jpi, jpj ) 594 575 CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj ) 595 576 ! 596 577 CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj ) 597 578 CALL tab_1d_2d( nbpb, hfx_spr , npb, hfx_spr_1d(1:nbpb) , jpi, jpj ) … … 612 593 ! 613 594 END SELECT 614 595 ! 615 596 END SUBROUTINE lim_thd_1d2d 616 597 … … 629 610 !!------------------------------------------------------------------- 630 611 INTEGER :: ios ! Local integer output status for namelist read 631 NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, & 632 & rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 612 !! 613 NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, & 614 & rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 633 615 & nn_monocat, ln_it_qnsice 634 616 !!------------------------------------------------------------------- -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r5517 r6140 41 41 CONTAINS 42 42 43 #if defined key_dimgout44 # include "limwri_dimg.h90"45 #else46 43 47 44 SUBROUTINE lim_wri( kindic ) … … 298 295 299 296 END SUBROUTINE lim_wri 300 #endif301 297 302 298
Note: See TracChangeset
for help on using the changeset viewer.