Changeset 6853 for branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
- Timestamp:
- 2016-08-08T11:02:18+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
r6746 r6853 31 31 USE iom 32 32 33 !!!clem 34 !! USE diawri 35 !!! 36 33 37 IMPLICIT NONE 34 38 PRIVATE … … 84 88 REAL(wp) :: ztmelts, zdh 85 89 INTEGER :: i_hemis, i_fill, jl0 86 REAL(wp) :: z test_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv90 REAL(wp) :: zarg, zV, zconv, zdv 87 91 REAL(wp), POINTER, DIMENSION(:,:) :: zswitch ! ice indicator 88 92 REAL(wp), POINTER, DIMENSION(:,:) :: zht_i_ini, zat_i_ini, zvt_i_ini !data from namelist or nc file 89 93 REAL(wp), POINTER, DIMENSION(:,:) :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 90 REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini, zv_i_ini !data by cattegories to fill 91 !-------------------------------------------------------------------- 92 93 CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini, za_i_ini, zv_i_ini ) 94 REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini !data by cattegories to fill 95 INTEGER , POINTER, DIMENSION(:) :: itest 96 !-------------------------------------------------------------------- 97 98 CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini, za_i_ini ) 94 99 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 ) 95 100 CALL wrk_alloc( jpi, jpj, zswitch ) 101 Call wrk_alloc( 4, itest ) 96 102 97 103 IF(lwp) WRITE(numout,*) … … 124 130 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 125 131 DO ji = 1, jpi 126 IF( ( sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN127 zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)! no ice132 IF( ( sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst * tmask(ji,jj,1) ) THEN 133 zswitch(ji,jj) = 0._wp ! no ice 128 134 ELSE 129 135 zswitch(ji,jj) = 1._wp * tmask(ji,jj,1) ! ice … … 153 159 DO ji = 1, jpi 154 160 IF( ff(ji,jj) >= 0._wp ) THEN 155 zht_i_ini(ji,jj) = rn_hti_ini_n 156 zht_s_ini(ji,jj) = rn_hts_ini_n 157 zat_i_ini(ji,jj) = rn_ati_ini_n 158 zts_u_ini(ji,jj) = rn_tmi_ini_n 159 zsm_i_ini(ji,jj) = rn_smi_ini_n 160 ztm_i_ini(ji,jj) = rn_tmi_ini_n 161 zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj) 162 zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj) 163 zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj) 164 zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj) 165 zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj) 166 ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj) 161 167 ELSE 162 zht_i_ini(ji,jj) = rn_hti_ini_s 163 zht_s_ini(ji,jj) = rn_hts_ini_s 164 zat_i_ini(ji,jj) = rn_ati_ini_s 165 zts_u_ini(ji,jj) = rn_tmi_ini_s 166 zsm_i_ini(ji,jj) = rn_smi_ini_s 167 ztm_i_ini(ji,jj) = rn_tmi_ini_s 168 zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj) 169 zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj) 170 zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj) 171 zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj) 172 zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj) 173 ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj) 168 174 ENDIF 169 175 END DO … … 181 187 zh_i_ini(:,:,:) = 0._wp 182 188 za_i_ini(:,:,:) = 0._wp 183 zv_i_ini(:,:,:) = 0._wp184 189 185 190 DO jj = 1, jpj … … 188 193 IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 189 194 190 ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 191 ! ztests = 0 192 193 DO i_fill = jpl, 1, -1 194 195 ! IF( ztests .NE. 4 ) THEN 196 IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 197 !---------------------------- 198 ! fill the i_fill categories 199 !---------------------------- 200 ! *** 1 category to fill 201 IF ( i_fill .EQ. 1 ) THEN 202 zh_i_ini(ji,jj, 1) = zht_i_ini(ji,jj) 203 za_i_ini(ji,jj, 1) = zat_i_ini(ji,jj) 204 zh_i_ini(ji,jj,2:jpl) = 0._wp 205 za_i_ini(ji,jj,2:jpl) = 0._wp 206 ELSE 207 208 ! *** >1 categores to fill 209 !--- Ice thicknesses in the i_fill - 1 first categories 210 DO jl = 1, i_fill - 1 211 zh_i_ini(ji,jj,jl) = hi_mean(jl) 212 END DO 195 !--- jl0: most likely index where cc will be maximum 196 jl0 = jpl 197 DO jl = 1, jpl 198 IF ( ( zht_i_ini(ji,jj) > hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 199 jl0 = jl 200 CYCLE 201 ENDIF 202 END DO 203 204 ! initialisation of tests 205 itest(:) = 0 206 207 i_fill = jpl + 1 !==================================== 208 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 209 ! iteration !==================================== 210 i_fill = i_fill - 1 211 212 ! initialisation of ice variables for each try 213 zh_i_ini(ji,jj,:) = 0._wp 214 za_i_ini(ji,jj,:) = 0._wp 215 itest(:) = 0 216 217 ! *** case very thin ice: fill only category 1 218 IF ( i_fill == 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 222 ! *** case ice is thicker: fill categories >1 223 ELSE 224 225 ! Fill ice thicknesses in the (i_fill-1) cat by hmean 226 DO jl = 1, i_fill-1 227 zh_i_ini(ji,jj,jl) = hi_mean(jl) 228 END DO 213 229 214 !--- jl0: most likely index where cc will be maximum 215 DO jl = 1, jpl 216 IF ( ( zht_i_ini(ji,jj) > hi_max(jl-1) ) .AND. & 217 & ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 218 jl0 = jl 219 ENDIF 220 END DO 221 jl0 = MIN(jl0, i_fill) 222 223 !--- Concentrations 224 za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 225 DO jl = 1, i_fill - 1 226 IF( jl .NE. jl0 )THEN 227 zsigma = 0.5 * zht_i_ini(ji,jj) 228 zarg = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / zsigma 229 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 230 ENDIF 231 END DO 232 233 zA = 0. ! sum of the areas in the jpl categories 234 DO jl = 1, i_fill - 1 235 zA = zA + za_i_ini(ji,jj,jl) 236 END DO 237 za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - zA ! ice conc in the last category 238 IF ( i_fill .LT. jpl ) za_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 239 240 !--- Ice thickness in the last category 241 zV = 0. ! sum of the volumes of the N-1 categories 242 DO jl = 1, i_fill - 1 243 zV = zV + za_i_ini(ji,jj,jl)*zh_i_ini(ji,jj,jl) 244 END DO 245 zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / za_i_ini(ji,jj,i_fill) 246 IF ( i_fill .LT. jpl ) zh_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 247 248 !--- volumes 249 zv_i_ini(ji,jj,:) = za_i_ini(ji,jj,:) * zh_i_ini(ji,jj,:) 250 IF ( i_fill .LT. jpl ) zv_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 251 252 ENDIF ! i_fill 253 254 !--------------------- 255 ! Compatibility tests 256 !--------------------- 257 ! Test 1: area conservation 258 zA_cons = SUM(za_i_ini(ji,jj,:)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons ) 259 IF ( zconv .LT. 1.0e-6 ) THEN 260 ztest_1 = 1 261 ELSE 262 ztest_1 = 0 263 ENDIF 264 265 ! Test 2: volume conservation 266 zV_cons = SUM(zv_i_ini(ji,jj,:)) 267 zconv = ABS(zvt_i_ini(ji,jj) - zV_cons) 268 269 IF( zconv .LT. 1.0e-6 ) THEN 270 ztest_2 = 1 271 ELSE 272 ztest_2 = 0 273 ENDIF 274 275 ! Test 3: thickness of the last category is in-bounds ? 276 IF ( zh_i_ini(ji,jj,i_fill) > hi_max(i_fill-1) ) THEN 277 ztest_3 = 1 278 ELSE 279 ztest_3 = 0 280 ENDIF 281 282 ! Test 4: positivity of ice concentrations 283 ztest_4 = 1 284 DO jl = 1, jpl 285 IF ( za_i_ini(ji,jj,jl) .LT. 0._wp ) THEN 286 ztest_4 = 0 230 !--- Concentrations 231 za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 232 DO jl = 1, i_fill - 1 233 IF( jl /= jl0 )THEN 234 zarg = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) ) 235 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 287 236 ENDIF 288 237 END DO 289 290 ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 291 292 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 293 294 END DO ! i_fill 295 296 IF(lwp) THEN 297 WRITE(numout,*) ' ztests : ', ztests 298 IF( ztests .NE. 4 )THEN 299 WRITE(numout,*) 300 WRITE(numout,*) ' !!!! ALERT !!! ' 301 WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 302 WRITE(numout,*) 303 WRITE(numout,*) ' *** ztests is not equal to 4 ' 304 WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 305 WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 306 WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 307 ENDIF ! ztests .NE. 4 238 239 ! Concentration in the last (i_fill) category 240 za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) ) 241 242 ! Ice thickness in the last (i_fill) category 243 zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) ) 244 zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 245 246 ! clem: correction if concentration of upper cat is greater than lower cat 247 ! (it should be a gaussian around jl0 but sometimes it is not) 248 IF ( jl0 /= jpl ) THEN 249 DO jl = jpl, jl0+1, -1 250 IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN 251 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl) 252 zh_i_ini(ji,jj,jl ) = 0._wp 253 za_i_ini(ji,jj,jl ) = 0._wp 254 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1) & 255 & + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 ) 256 END IF 257 ENDDO 258 ENDIF 259 260 ENDIF ! case ice is thick or thin 261 262 !--------------------- 263 ! Compatibility tests 264 !--------------------- 265 ! Test 1: area conservation 266 zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) ) 267 IF ( zconv < epsi06 ) itest(1) = 1 268 269 ! Test 2: volume conservation 270 zconv = ABS( zat_i_ini(ji,jj) * zht_i_ini(ji,jj) & 271 & - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) ) 272 IF ( zconv < epsi06 ) itest(2) = 1 273 274 ! Test 3: thickness of the last category is in-bounds ? 275 IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 276 277 ! Test 4: positivity of ice concentrations 278 itest(4) = 1 279 DO jl = 1, i_fill 280 IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0 281 END DO 282 ! !============================ 283 END DO ! end iteration on categories 284 ! !============================ 285 286 IF( lwp .AND. SUM(itest) /= 4 ) THEN 287 WRITE(numout,*) 288 WRITE(numout,*) ' !!!! ALERT itest is not equal to 4 !!! ' 289 WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 290 WRITE(numout,*) 291 WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:) 292 WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 293 WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 308 294 ENDIF 309 310 ENDIF ! zat_i_ini(ji,jj) > 0._wp .AND. zh m_i_ini(ji,jj) > 0._wp311 312 ENDDO 313 ENDDO 295 296 ENDIF ! zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp 297 298 ENDDO 299 ENDDO 314 300 315 301 !--------------------------------------------------------------------- … … 468 454 sxyage (:,:,:) = 0._wp 469 455 470 471 CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini, za_i_ini, zv_i_ini ) 456 !!!clem 457 !! ! Output the initial state and forcings 458 !! CALL dia_wri_state( 'output.init', nit000 ) 459 !!! 460 461 CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini, za_i_ini ) 472 462 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 ) 473 463 CALL wrk_dealloc( jpi, jpj, zswitch ) 464 Call wrk_dealloc( 4, itest ) 474 465 475 466 END SUBROUTINE lim_istate
Note: See TracChangeset
for help on using the changeset viewer.
