- Timestamp:
- 2020-12-03T12:20:38+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icectl.F90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette @13292sette10 ^/utils/CI/sette_wave@13990 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icectl.F90
r13295 r14037 43 43 PUBLIC ice_prt 44 44 PUBLIC ice_prt3D 45 PUBLIC ice_drift_wri 46 PUBLIC ice_drift_init 45 47 46 48 ! thresold rates for conservation … … 49 51 REAL(wp), PARAMETER :: zchk_s = 2.5e-6 ! g/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg) 50 52 REAL(wp), PARAMETER :: zchk_t = 7.5e-2 ! W/m2 <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg) 53 54 ! for drift outputs 55 CHARACTER(LEN=50) :: clname="icedrift_diagnostics.ascii" ! ascii filename 56 INTEGER :: numicedrift ! outfile unit 57 REAL(wp) :: rdiag_icemass, rdiag_icesalt, rdiag_iceheat 58 REAL(wp) :: rdiag_adv_icemass, rdiag_adv_icesalt, rdiag_adv_iceheat 51 59 52 60 !! * Substitutions … … 77 85 !! 78 86 REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, & 79 & zdiag_vmin, zdiag_amin, zdiag_amax, zdiag_eimin, zdiag_esmin, zdiag_smin 87 & zdiag_vimin, zdiag_vsmin, zdiag_vpmin, zdiag_vlmin, zdiag_aimin, zdiag_aimax, & 88 & zdiag_eimin, zdiag_esmin, zdiag_simin 80 89 REAL(wp) :: zvtrp, zetrp 81 90 REAL(wp) :: zarea … … 84 93 IF( icount == 0 ) THEN 85 94 86 pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos , dim=3 ) * e1e2t )95 pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) 87 96 pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) 88 97 pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ) … … 104 113 105 114 ! -- mass diag -- ! 106 zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_Dt_ice & 115 zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) & 116 & - pdiag_v ) * r1_Dt_ice & 107 117 & + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + & 108 118 & wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + & … … 124 134 125 135 ! -- min/max diag -- ! 126 zdiag_amax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 127 zdiag_vmin = glob_min( 'icectl', v_i ) 128 zdiag_amin = glob_min( 'icectl', a_i ) 129 zdiag_smin = glob_min( 'icectl', sv_i ) 136 zdiag_aimax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 137 zdiag_vimin = glob_min( 'icectl', v_i ) 138 zdiag_vsmin = glob_min( 'icectl', v_s ) 139 zdiag_vpmin = glob_min( 'icectl', v_ip ) 140 zdiag_vlmin = glob_min( 'icectl', v_il ) 141 zdiag_aimin = glob_min( 'icectl', a_i ) 142 zdiag_simin = glob_min( 'icectl', sv_i ) 130 143 zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 131 144 zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 132 145 133 146 ! -- advection scheme is conservative? -- ! 134 zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 (only for Prather)135 zetrp = glob_sum( 'icectl', ( diag_trp_ei + diag_trp_es ) * e1e2t ) ! must be close to 0 (only for Prather)147 zvtrp = glob_sum( 'icectl', diag_adv_mass * e1e2t ) 148 zetrp = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 136 149 137 150 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 147 160 & WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 148 161 ! check negative values 149 IF( zdiag_vmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_i < 0 = ',zdiag_vmin 150 IF( zdiag_amin < 0. ) WRITE(numout,*) cd_routine,' : violation a_i < 0 = ',zdiag_amin 151 IF( zdiag_smin < 0. ) WRITE(numout,*) cd_routine,' : violation s_i < 0 = ',zdiag_smin 152 IF( zdiag_eimin < 0. ) WRITE(numout,*) cd_routine,' : violation e_i < 0 = ',zdiag_eimin 153 IF( zdiag_esmin < 0. ) WRITE(numout,*) cd_routine,' : violation e_s < 0 = ',zdiag_esmin 162 IF( zdiag_vimin < 0. ) WRITE(numout,*) cd_routine,' : violation v_i < 0 = ',zdiag_vimin 163 IF( zdiag_vsmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_s < 0 = ',zdiag_vsmin 164 IF( zdiag_vpmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_ip < 0 = ',zdiag_vpmin 165 IF( zdiag_vlmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_il < 0 = ',zdiag_vlmin 166 IF( zdiag_aimin < 0. ) WRITE(numout,*) cd_routine,' : violation a_i < 0 = ',zdiag_aimin 167 IF( zdiag_simin < 0. ) WRITE(numout,*) cd_routine,' : violation s_i < 0 = ',zdiag_simin 168 IF( zdiag_eimin < 0. ) WRITE(numout,*) cd_routine,' : violation e_i < 0 = ',zdiag_eimin 169 IF( zdiag_esmin < 0. ) WRITE(numout,*) cd_routine,' : violation e_s < 0 = ',zdiag_esmin 154 170 ! check maximum ice concentration 155 IF( zdiag_a max >MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) &156 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_a max171 IF( zdiag_aimax>MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 172 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_aimax 157 173 ! check if advection scheme is conservative 158 ! only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc) 159 ! so the formulation for conservation is different (and not coded) 160 ! it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now 161 !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 162 ! & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 174 IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 175 & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 176 IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 177 & WRITE(numout,*) cd_routine,' : violation adv scheme [J] = ',zetrp * rDt_ice 163 178 ENDIF 164 179 ! … … 186 201 ! water flux 187 202 ! -- mass diag -- ! 188 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) 203 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + wfx_pnd & 204 & + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ) 189 205 190 206 ! -- salt diag -- ! 191 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t )207 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) 192 208 193 209 ! -- heat diag -- ! 194 ! clem: not the good formulation 195 !!zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr & 196 !! & ) * e1e2t ) 210 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 211 ! equivalent to this: 212 !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 213 !! & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr & 214 !! & ) * e1e2t ) 197 215 198 216 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 204 222 IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 205 223 & WRITE(numout,*) cd_routine,' : violation salt cons. [g] = ',zdiag_salt * rDt_ice 206 !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 224 IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 225 & WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 207 226 ENDIF 208 227 ! … … 234 253 IF( icount == 0 ) THEN 235 254 236 pdiag_v = SUM( v_i * rhoi + v_s * rhos , dim=3 )237 pdiag_s = SUM( sv_i * rhoi 255 pdiag_v = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) 256 pdiag_s = SUM( sv_i * rhoi , dim=3 ) 238 257 pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) 239 258 … … 250 269 251 270 ! -- mass diag -- ! 252 zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos , dim=3 ) - pdiag_v ) * r1_Dt_ice&271 zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) - pdiag_v ) * r1_Dt_ice & 253 272 & + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 254 273 & wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) & … … 341 360 CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 ) ! 342 361 CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 ) ! 362 ! mean state 363 CALL iom_rstput( 0, 0, inum, 'icecon' , SUM(a_i ,dim=3) , ktype = jp_r8 ) ! 364 CALL iom_rstput( 0, 0, inum, 'icevol' , SUM(v_i ,dim=3) , ktype = jp_r8 ) ! 365 CALL iom_rstput( 0, 0, inum, 'snwvol' , SUM(v_s ,dim=3) , ktype = jp_r8 ) ! 366 CALL iom_rstput( 0, 0, inum, 'pndvol' , SUM(v_ip,dim=3) , ktype = jp_r8 ) ! 367 CALL iom_rstput( 0, 0, inum, 'lidvol' , SUM(v_il,dim=3) , ktype = jp_r8 ) ! 343 368 344 369 CALL iom_close( inum ) … … 350 375 !! *** ROUTINE ice_ctl *** 351 376 !! 352 !! ** Purpose : Alerts in case of model crash377 !! ** Purpose : control checks 353 378 !!------------------------------------------------------------------- 354 379 INTEGER, INTENT(in) :: kt ! ocean time step 355 INTEGER :: ji, jj, jk, jl ! dummy loop indices 356 INTEGER :: inb_altests ! number of alert tests (max 20) 357 INTEGER :: ialert_id ! number of the current alert 358 REAL(wp) :: ztmelts ! ice layer melting point 380 INTEGER :: ja, ji, jj, jk, jl ! dummy loop indices 381 INTEGER :: ialert_id ! number of the current alert 382 REAL(wp) :: ztmelts ! ice layer melting point 359 383 CHARACTER (len=30), DIMENSION(20) :: cl_alname ! name of alert 360 384 INTEGER , DIMENSION(20) :: inb_alp ! number of alerts positive 361 385 !!------------------------------------------------------------------- 362 363 inb_altests = 10 364 inb_alp(:) = 0 365 366 ! Alert if incompatible volume and concentration 367 ialert_id = 2 ! reference number of this alert 368 cl_alname(ialert_id) = ' Incompat vol and con ' ! name of the alert 386 inb_alp(:) = 0 387 ialert_id = 0 388 389 ! Alert if very high salinity 390 ialert_id = ialert_id + 1 ! reference number of this alert 391 cl_alname(ialert_id) = ' Very high salinity ' ! name of the alert 369 392 DO jl = 1, jpl 370 393 DO_2D( 1, 1, 1, 1 ) 371 IF( v_i(ji,jj,jl) /= 0._wp .AND. a_i(ji,jj,jl) == 0._wp ) THEN 372 WRITE(numout,*) ' ALERTE 2 : Incompatible volume and concentration ' 373 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 394 IF( v_i(ji,jj,jl) > epsi10 ) THEN 395 IF( sv_i(ji,jj,jl) / v_i(ji,jj,jl) > rn_simax ) THEN 396 WRITE(numout,*) ' ALERTE : Very high salinity ',sv_i(ji,jj,jl)/v_i(ji,jj,jl) 397 WRITE(numout,*) ' at i,j,l = ',ji,jj,jl 398 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 399 ENDIF 374 400 ENDIF 375 401 END_2D 376 402 END DO 377 403 378 ! Alerte if very thick ice 379 ialert_id = 3 ! reference number of this alert 380 cl_alname(ialert_id) = ' Very thick ice ' ! name of the alert 381 jl = jpl 382 DO_2D( 1, 1, 1, 1 ) 383 IF( h_i(ji,jj,jl) > 50._wp ) THEN 384 WRITE(numout,*) ' ALERTE 3 : Very thick ice' 385 !CALL ice_prt( kt, ji, jj, 2, ' ALERTE 3 : Very thick ice ' ) 386 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 387 ENDIF 388 END_2D 389 390 ! Alert if very fast ice 391 ialert_id = 4 ! reference number of this alert 392 cl_alname(ialert_id) = ' Very fast ice ' ! name of the alert 393 DO_2D( 1, 1, 1, 1 ) 394 IF( MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 2. .AND. & 395 & at_i(ji,jj) > 0._wp ) THEN 396 WRITE(numout,*) ' ALERTE 4 : Very fast ice' 397 !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 4 : Very fast ice ' ) 398 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 399 ENDIF 400 END_2D 401 402 ! Alert on salt flux 403 ialert_id = 5 ! reference number of this alert 404 cl_alname(ialert_id) = ' High salt flux ' ! name of the alert 405 DO_2D( 1, 1, 1, 1 ) 406 IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN ! = 1 psu/day for 1m ocean depth 407 WRITE(numout,*) ' ALERTE 5 : High salt flux' 408 !CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 : High salt flux ' ) 409 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 410 ENDIF 411 END_2D 412 413 ! Alert if there is ice on continents 414 ialert_id = 6 ! reference number of this alert 415 cl_alname(ialert_id) = ' Ice on continents ' ! name of the alert 416 DO_2D( 1, 1, 1, 1 ) 417 IF( tmask(ji,jj,1) <= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 418 WRITE(numout,*) ' ALERTE 6 : Ice on continents' 419 !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 6 : Ice on continents ' ) 420 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 421 ENDIF 422 END_2D 423 424 ! 425 ! ! Alert if very fresh ice 426 ialert_id = 7 ! reference number of this alert 427 cl_alname(ialert_id) = ' Very fresh ice ' ! name of the alert 404 ! Alert if very low salinity 405 ialert_id = ialert_id + 1 ! reference number of this alert 406 cl_alname(ialert_id) = ' Very low salinity ' ! name of the alert 428 407 DO jl = 1, jpl 429 408 DO_2D( 1, 1, 1, 1 ) 430 IF( s_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 431 WRITE(numout,*) ' ALERTE 7 : Very fresh ice' 432 ! CALL ice_prt(kt,ji,jj,1, ' ALERTE 7 : Very fresh ice ' ) 433 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 409 IF( v_i(ji,jj,jl) > epsi10 ) THEN 410 IF( sv_i(ji,jj,jl) / v_i(ji,jj,jl) < rn_simin ) THEN 411 WRITE(numout,*) ' ALERTE : Very low salinity ',sv_i(ji,jj,jl),v_i(ji,jj,jl) 412 WRITE(numout,*) ' at i,j,l = ',ji,jj,jl 413 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 414 ENDIF 434 415 ENDIF 435 416 END_2D 436 417 END DO 437 ! 438 ! Alert if qns very big 439 ialert_id = 8 ! reference number of this alert 440 cl_alname(ialert_id) = ' fnsolar very big ' ! name of the alert 441 DO_2D( 1, 1, 1, 1 ) 442 IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 443 ! 444 WRITE(numout,*) ' ALERTE 8 : Very high non-solar heat flux' 445 !CALL ice_prt( kt, ji, jj, 2, ' ') 446 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 447 ! 448 ENDIF 449 END_2D 450 !+++++ 451 452 ! ! Alert if too old ice 453 ialert_id = 9 ! reference number of this alert 454 cl_alname(ialert_id) = ' Very old ice ' ! name of the alert 455 DO jl = 1, jpl 456 DO_2D( 1, 1, 1, 1 ) 457 IF ( ( ( ABS( o_i(ji,jj,jl) ) > rDt_ice ) .OR. & 458 ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 459 ( a_i(ji,jj,jl) > 0._wp ) ) THEN 460 WRITE(numout,*) ' ALERTE 9 : Wrong ice age' 461 !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 9 : Wrong ice age ') 462 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 463 ENDIF 464 END_2D 465 END DO 466 467 ! Alert if very warm ice 468 ialert_id = 10 ! reference number of this alert 469 cl_alname(ialert_id) = ' Very warm ice ' ! name of the alert 470 inb_alp(ialert_id) = 0 418 419 ! Alert if very cold ice 420 ialert_id = ialert_id + 1 ! reference number of this alert 421 cl_alname(ialert_id) = ' Very cold ice ' ! name of the alert 471 422 DO jl = 1, jpl 472 423 DO_3D( 1, 1, 1, 1, 1, nlay_i ) 473 424 ztmelts = -rTmlt * sz_i(ji,jj,jk,jl) + rt0 474 IF( t_i(ji,jj,jk,jl) > ztmelts .AND. v_i(ji,jj,jl) > 1.e-10 &475 & .AND. a_i(ji,jj,jl) > 0._wp ) THEN476 WRITE(numout,*) ' ALERTE 10 : Very warm ice'425 IF( t_i(ji,jj,jk,jl) < -50.+rt0 .AND. v_i(ji,jj,jl) > epsi10 ) THEN 426 WRITE(numout,*) ' ALERTE : Very cold ice ',(t_i(ji,jj,jk,jl)-rt0) 427 WRITE(numout,*) ' at i,j,k,l = ',ji,jj,jk,jl 477 428 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 478 429 ENDIF 479 430 END_3D 480 431 END DO 432 433 ! Alert if very warm ice 434 ialert_id = ialert_id + 1 ! reference number of this alert 435 cl_alname(ialert_id) = ' Very warm ice ' ! name of the alert 436 DO jl = 1, jpl 437 DO_3D( 1, 1, 1, 1, 1, nlay_i ) 438 ztmelts = -rTmlt * sz_i(ji,jj,jk,jl) + rt0 439 IF( t_i(ji,jj,jk,jl) > ztmelts .AND. v_i(ji,jj,jl) > epsi10 ) THEN 440 WRITE(numout,*) ' ALERTE : Very warm ice',(t_i(ji,jj,jk,jl)-rt0) 441 WRITE(numout,*) ' at i,j,k,l = ',ji,jj,jk,jl 442 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 443 ENDIF 444 END_3D 445 END DO 446 447 ! Alerte if very thick ice 448 ialert_id = ialert_id + 1 ! reference number of this alert 449 cl_alname(ialert_id) = ' Very thick ice ' ! name of the alert 450 jl = jpl 451 DO_2D( 1, 1, 1, 1 ) 452 IF( h_i(ji,jj,jl) > 50._wp ) THEN 453 WRITE(numout,*) ' ALERTE : Very thick ice ',h_i(ji,jj,jl) 454 WRITE(numout,*) ' at i,j,l = ',ji,jj,jl 455 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 456 ENDIF 457 END_2D 458 459 ! Alerte if very thin ice 460 ialert_id = ialert_id + 1 ! reference number of this alert 461 cl_alname(ialert_id) = ' Very thin ice ' ! name of the alert 462 jl = 1 463 DO_2D( 1, 1, 1, 1 ) 464 IF( h_i(ji,jj,jl) < rn_himin ) THEN 465 WRITE(numout,*) ' ALERTE : Very thin ice ',h_i(ji,jj,jl) 466 WRITE(numout,*) ' at i,j,l = ',ji,jj,jl 467 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 468 ENDIF 469 END_2D 470 471 ! Alert if very fast ice 472 ialert_id = ialert_id + 1 ! reference number of this alert 473 cl_alname(ialert_id) = ' Very fast ice ' ! name of the alert 474 DO_2D( 1, 1, 1, 1 ) 475 IF( MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 2. ) THEN 476 WRITE(numout,*) ' ALERTE : Very fast ice ',MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) 477 WRITE(numout,*) ' at i,j = ',ji,jj 478 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 479 ENDIF 480 END_2D 481 482 ! Alert if there is ice on continents 483 ialert_id = ialert_id + 1 ! reference number of this alert 484 cl_alname(ialert_id) = ' Ice on continents ' ! name of the alert 485 DO_2D( 1, 1, 1, 1 ) 486 IF( tmask(ji,jj,1) == 0._wp .AND. ( at_i(ji,jj) > 0._wp .OR. vt_i(ji,jj) > 0._wp ) ) THEN 487 WRITE(numout,*) ' ALERTE : Ice on continents ',at_i(ji,jj),vt_i(ji,jj) 488 WRITE(numout,*) ' at i,j = ',ji,jj 489 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 490 ENDIF 491 END_2D 492 493 ! Alert if incompatible ice concentration and volume 494 ialert_id = ialert_id + 1 ! reference number of this alert 495 cl_alname(ialert_id) = ' Incompatible ice conc and vol ' ! name of the alert 496 DO_2D( 1, 1, 1, 1 ) 497 IF( ( vt_i(ji,jj) == 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. & 498 & ( vt_i(ji,jj) > 0._wp .AND. at_i(ji,jj) == 0._wp ) ) THEN 499 WRITE(numout,*) ' ALERTE : Incompatible ice conc and vol ',at_i(ji,jj),vt_i(ji,jj) 500 WRITE(numout,*) ' at i,j = ',ji,jj 501 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 502 ENDIF 503 END_2D 481 504 482 505 ! sum of the alerts on all processors 483 506 IF( lk_mpp ) THEN 484 DO ialert_id = 1, inb_altests485 CALL mpp_sum('icectl', inb_alp( ialert_id))507 DO ja = 1, ialert_id 508 CALL mpp_sum('icectl', inb_alp(ja)) 486 509 END DO 487 510 ENDIF … … 489 512 ! print alerts 490 513 IF( lwp ) THEN 491 ialert_id = 1 ! reference number of this alert492 cl_alname(ialert_id) = ' NO alerte 1 ' ! name of the alert493 514 WRITE(numout,*) ' time step ',kt 494 515 WRITE(numout,*) ' All alerts at the end of ice model ' 495 DO ialert_id = 1, inb_altests496 WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! '516 DO ja = 1, ialert_id 517 WRITE(numout,*) ja, cl_alname(ja)//' : ', inb_alp(ja), ' times ! ' 497 518 END DO 498 519 ENDIF … … 543 564 WRITE(numout,*) ' v_ice(i ,j) : ', v_ice(ji,jj) 544 565 WRITE(numout,*) ' strength : ', strength(ji,jj) 545 WRITE(numout,*)546 566 WRITE(numout,*) ' - Cell values ' 547 567 WRITE(numout,*) ' ~~~~~~~~~~~ ' … … 552 572 DO jl = 1, jpl 553 573 WRITE(numout,*) ' - Category (', jl,')' 574 WRITE(numout,*) ' ~~~~~~~~~~~ ' 554 575 WRITE(numout,*) ' a_i : ', a_i(ji,jj,jl) 555 576 WRITE(numout,*) ' h_i : ', h_i(ji,jj,jl) … … 588 609 WRITE(numout,*) ' v_ice(i ,j) : ', v_ice(ji,jj) 589 610 WRITE(numout,*) ' strength : ', strength(ji,jj) 590 WRITE(numout,*) ' u_ice_b : ', u_ice_b(ji,jj) , ' v_ice_b : ', v_ice_b(ji,jj)591 611 WRITE(numout,*) 592 612 … … 605 625 WRITE(numout,*) ' e_snow : ', e_s(ji,jj,1,jl) , ' e_snow_b : ', e_s_b(ji,jj,1,jl) 606 626 WRITE(numout,*) ' sv_i : ', sv_i(ji,jj,jl) , ' sv_i_b : ', sv_i_b(ji,jj,jl) 607 WRITE(numout,*) ' oa_i : ', oa_i(ji,jj,jl) , ' oa_i_b : ', oa_i_b(ji,jj,jl)608 627 END DO !jl 609 628 … … 713 732 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' v_i : ') 714 733 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' v_s : ') 715 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl) , clinfo1= ' e_i1 : ')716 734 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' e_snow : ') 717 735 CALL prt_ctl(tab2d_1=sv_i (:,:,jl) , clinfo1= ' sv_i : ') … … 721 739 CALL prt_ctl_info(' - Layer : ', ivar=jk) 722 740 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' t_i : ') 741 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' e_i : ') 723 742 END DO 724 743 END DO … … 731 750 732 751 END SUBROUTINE ice_prt3D 752 753 754 SUBROUTINE ice_drift_wri( kt ) 755 !!------------------------------------------------------------------- 756 !! *** ROUTINE ice_drift_wri *** 757 !! 758 !! ** Purpose : conservation of mass, salt and heat 759 !! write the drift in a ascii file at each time step 760 !! and the total run drifts 761 !!------------------------------------------------------------------- 762 INTEGER, INTENT(in) :: kt ! ice time-step index 763 ! 764 INTEGER :: ji, jj 765 REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 766 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_mass2D, zdiag_salt2D, zdiag_heat2D 767 !!------------------------------------------------------------------- 768 ! 769 IF( kt == nit000 .AND. lwp ) THEN 770 WRITE(numout,*) 771 WRITE(numout,*) 'ice_drift_wri: sea-ice drifts' 772 WRITE(numout,*) '~~~~~~~~~~~~~' 773 ENDIF 774 ! 775 ! 2D budgets (must be close to 0) 776 IF( iom_use('icedrift_mass') .OR. iom_use('icedrift_salt') .OR. iom_use('icedrift_heat') ) THEN 777 DO_2D( 1, 1, 1, 1 ) 778 zdiag_mass2D(ji,jj) = wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_spr(ji,jj) + wfx_sub(ji,jj) & 779 & + diag_vice(ji,jj) + diag_vsnw(ji,jj) - diag_adv_mass(ji,jj) 780 zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) 781 zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj) 782 END_2D 783 ! 784 ! write outputs 785 CALL iom_put( 'icedrift_mass', zdiag_mass2D ) 786 CALL iom_put( 'icedrift_salt', zdiag_salt2D ) 787 CALL iom_put( 'icedrift_heat', zdiag_heat2D ) 788 ENDIF 789 790 ! -- mass diag -- ! 791 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub & 792 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rDt_ice 793 zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 794 795 ! -- salt diag -- ! 796 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rDt_ice * 1.e-3 797 zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 798 799 ! -- heat diag -- ! 800 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 801 zdiag_adv_heat = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 802 803 ! ! write out to file 804 IF( lwp ) THEN 805 ! check global drift (must be close to 0) 806 WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift [kg]', zdiag_mass 807 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift [kg]', zdiag_salt 808 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift [W] ', zdiag_heat 809 ! check drift from advection scheme (can be /=0 with bdy but not sure why) 810 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'mass drift adv [kg]', zdiag_adv_mass 811 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift adv [kg]', zdiag_adv_salt 812 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift adv [W] ', zdiag_adv_heat 813 ENDIF 814 ! ! drifts 815 rdiag_icemass = rdiag_icemass + zdiag_mass 816 rdiag_icesalt = rdiag_icesalt + zdiag_salt 817 rdiag_iceheat = rdiag_iceheat + zdiag_heat 818 rdiag_adv_icemass = rdiag_adv_icemass + zdiag_adv_mass 819 rdiag_adv_icesalt = rdiag_adv_icesalt + zdiag_adv_salt 820 rdiag_adv_iceheat = rdiag_adv_iceheat + zdiag_adv_heat 821 ! 822 ! ! output drifts and close ascii file 823 IF( kt == nitend - nn_fsbc + 1 .AND. lwp ) THEN 824 ! to ascii file 825 WRITE(numicedrift,*) '******************************************' 826 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift [kg]', rdiag_icemass 827 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift adv [kg]', rdiag_adv_icemass 828 WRITE(numicedrift,*) '******************************************' 829 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift [kg]', rdiag_icesalt 830 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift adv [kg]', rdiag_adv_icesalt 831 WRITE(numicedrift,*) '******************************************' 832 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift [W] ', rdiag_iceheat 833 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift adv [W] ', rdiag_adv_iceheat 834 CLOSE( numicedrift ) 835 ! 836 ! to ocean output 837 WRITE(numout,*) 838 WRITE(numout,*) 'ice_drift_wri: ice drifts information for the run ' 839 WRITE(numout,*) '~~~~~~~~~~~~~' 840 ! check global drift (must be close to 0) 841 WRITE(numout,*) ' sea-ice mass drift [kg] = ', rdiag_icemass 842 WRITE(numout,*) ' sea-ice salt drift [kg] = ', rdiag_icesalt 843 WRITE(numout,*) ' sea-ice heat drift [W] = ', rdiag_iceheat 844 ! check drift from advection scheme (can be /=0 with bdy but not sure why) 845 WRITE(numout,*) ' sea-ice mass drift adv [kg] = ', rdiag_adv_icemass 846 WRITE(numout,*) ' sea-ice salt drift adv [kg] = ', rdiag_adv_icesalt 847 WRITE(numout,*) ' sea-ice heat drift adv [W] = ', rdiag_adv_iceheat 848 ENDIF 849 ! 850 END SUBROUTINE ice_drift_wri 851 852 SUBROUTINE ice_drift_init 853 !!---------------------------------------------------------------------- 854 !! *** ROUTINE ice_drift_init *** 855 !! 856 !! ** Purpose : create output file, initialise arrays 857 !!---------------------------------------------------------------------- 858 ! 859 IF( .NOT.ln_icediachk ) RETURN ! exit 860 ! 861 IF(lwp) THEN 862 WRITE(numout,*) 863 WRITE(numout,*) 'ice_drift_init: Output ice drifts to ',TRIM(clname), ' file' 864 WRITE(numout,*) '~~~~~~~~~~~~~' 865 WRITE(numout,*) 866 ! 867 ! create output ascii file 868 CALL ctl_opn( numicedrift, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, narea ) 869 WRITE(numicedrift,*) 'Timestep Drifts' 870 WRITE(numicedrift,*) '******************************************' 871 ENDIF 872 ! 873 rdiag_icemass = 0._wp 874 rdiag_icesalt = 0._wp 875 rdiag_iceheat = 0._wp 876 rdiag_adv_icemass = 0._wp 877 rdiag_adv_icesalt = 0._wp 878 rdiag_adv_iceheat = 0._wp 879 ! 880 END SUBROUTINE ice_drift_init 733 881 734 882 #else
Note: See TracChangeset
for help on using the changeset viewer.
