Changeset 14072 for NEMO/trunk/src/ICE/icectl.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icectl.F90 (modified) (35 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icectl.F90
r14005 r14072 12 12 !! 'key_si3' SI3 sea-ice model 13 13 !!---------------------------------------------------------------------- 14 !! ice_cons_hsm : conservation tests on heat, salt and mass during a time step (global) 14 !! ice_cons_hsm : conservation tests on heat, salt and mass during a time step (global) 15 15 !! ice_cons_final : conservation tests on heat, salt and mass at end of time step (global) 16 16 !! ice_cons2D : conservation tests on heat, salt and mass at each gridcell … … 55 55 CHARACTER(LEN=50) :: clname="icedrift_diagnostics.ascii" ! ascii filename 56 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 59 57 REAL(wp) :: rdiag_icemass, rdiag_icesalt, rdiag_iceheat 58 REAL(wp) :: rdiag_adv_icemass, rdiag_adv_icesalt, rdiag_adv_iceheat 59 60 60 !! * Substitutions 61 61 # include "do_loop_substitute.h90" … … 77 77 !! It prints in ocean.output if there is a violation of conservation at each time-step 78 78 !! The thresholds (zchk_m, zchk_s, zchk_t) determine violations 79 !! For salt and heat thresholds, ice is considered to have a salinity of 10 80 !! and a heat content of 3e5 J/kg (=latent heat of fusion) 79 !! For salt and heat thresholds, ice is considered to have a salinity of 10 80 !! and a heat content of 3e5 J/kg (=latent heat of fusion) 81 81 !!------------------------------------------------------------------- 82 82 INTEGER , INTENT(in) :: icount ! called at: =0 the begining of the routine, =1 the end … … 148 148 zetrp = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 149 149 150 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 150 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 151 151 zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 152 152 … … 191 191 !! It prints in ocean.output if there is a violation of conservation at each time-step 192 192 !! The thresholds (zchk_m, zchk_s, zchk_t) determine the violations 193 !! For salt and heat thresholds, ice is considered to have a salinity of 10 194 !! and a heat content of 3e5 J/kg (=latent heat of fusion) 193 !! For salt and heat thresholds, ice is considered to have a salinity of 10 194 !! and a heat content of 3e5 J/kg (=latent heat of fusion) 195 195 !!------------------------------------------------------------------- 196 196 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine … … 214 214 !! & ) * e1e2t ) 215 215 216 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 216 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 217 217 zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 218 218 … … 243 243 !! 244 244 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_mass, zdiag_salt, zdiag_heat, & 245 & zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax 245 & zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax 246 246 INTEGER :: jl, jk 247 247 LOGICAL :: ll_stop_m = .FALSE. … … 261 261 & wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr 262 262 ! salt flux 263 pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam 263 pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam 264 264 ! heat flux 265 pdiag_ft = hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 265 pdiag_ft = hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 266 266 & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr 267 267 … … 283 283 ! -- heat diag -- ! 284 284 zdiag_heat = ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) - pdiag_t ) * r1_Dt_ice & 285 & + ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 285 & + ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 286 286 & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) & 287 287 & - pdiag_ft … … 324 324 IF( ll_stop_s ) CALL ctl_stop( 'STOP', cd_routine//': ice salt conservation issue' ) 325 325 IF( ll_stop_t ) CALL ctl_stop( 'STOP', cd_routine//': ice heat conservation issue' ) 326 326 327 327 ENDIF 328 328 … … 332 332 !!--------------------------------------------------------------------- 333 333 !! *** ROUTINE ice_cons_wri *** 334 !! 335 !! ** Purpose : create a NetCDF file named cdfile_name which contains 334 !! 335 !! ** Purpose : create a NetCDF file named cdfile_name which contains 336 336 !! the instantaneous fields when conservation issue occurs 337 337 !! … … 340 340 CHARACTER(len=*), INTENT( in ) :: cdfile_name ! name of the file created 341 341 REAL(wp), DIMENSION(:,:), INTENT( in ) :: pdiag_mass, pdiag_salt, pdiag_heat, & 342 & pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin !!, pdiag_amax 342 & pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin !!, pdiag_amax 343 343 !! 344 344 INTEGER :: inum 345 345 !!---------------------------------------------------------------------- 346 ! 346 ! 347 347 IF(lwp) WRITE(numout,*) 348 348 IF(lwp) WRITE(numout,*) 'ice_cons_wri : single instantaneous ice state' 349 349 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ named :', cdfile_name, '...nc' 350 IF(lwp) WRITE(numout,*) 350 IF(lwp) WRITE(numout,*) 351 351 352 352 CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' ) 353 353 354 354 CALL iom_rstput( 0, 0, inum, 'cons_mass', pdiag_mass(:,:) , ktype = jp_r8 ) ! ice mass spurious lost/gain 355 355 CALL iom_rstput( 0, 0, inum, 'cons_salt', pdiag_salt(:,:) , ktype = jp_r8 ) ! ice salt spurious lost/gain 356 356 CALL iom_rstput( 0, 0, inum, 'cons_heat', pdiag_heat(:,:) , ktype = jp_r8 ) ! ice heat spurious lost/gain 357 357 ! other diags 358 CALL iom_rstput( 0, 0, inum, 'aneg_count', pdiag_amin(:,:) , ktype = jp_r8 ) ! 359 CALL iom_rstput( 0, 0, inum, 'vneg_count', pdiag_vmin(:,:) , ktype = jp_r8 ) ! 360 CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 ) ! 361 CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 ) ! 358 CALL iom_rstput( 0, 0, inum, 'aneg_count', pdiag_amin(:,:) , ktype = jp_r8 ) ! 359 CALL iom_rstput( 0, 0, inum, 'vneg_count', pdiag_vmin(:,:) , ktype = jp_r8 ) ! 360 CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 ) ! 361 CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 ) ! 362 362 ! mean state 363 363 CALL iom_rstput( 0, 0, inum, 'icecon' , SUM(a_i ,dim=3) , ktype = jp_r8 ) ! … … 366 366 CALL iom_rstput( 0, 0, inum, 'pndvol' , SUM(v_ip,dim=3) , ktype = jp_r8 ) ! 367 367 CALL iom_rstput( 0, 0, inum, 'lidvol' , SUM(v_il,dim=3) , ktype = jp_r8 ) ! 368 368 369 369 CALL iom_close( inum ) 370 370 371 371 END SUBROUTINE ice_cons_wri 372 372 373 373 SUBROUTINE ice_ctl( kt ) 374 374 !!------------------------------------------------------------------- 375 !! *** ROUTINE ice_ctl *** 376 !! 375 !! *** ROUTINE ice_ctl *** 376 !! 377 377 !! ** Purpose : control checks 378 378 !!------------------------------------------------------------------- … … 386 386 inb_alp(:) = 0 387 387 ialert_id = 0 388 388 389 389 ! Alert if very high salinity 390 390 ialert_id = ialert_id + 1 ! reference number of this alert … … 430 430 END_3D 431 431 END DO 432 432 433 433 ! Alert if very warm ice 434 434 ialert_id = ialert_id + 1 ! reference number of this alert … … 444 444 END_3D 445 445 END DO 446 446 447 447 ! Alerte if very thick ice 448 448 ialert_id = ialert_id + 1 ! reference number of this alert 449 449 cl_alname(ialert_id) = ' Very thick ice ' ! name of the alert 450 jl = jpl 450 jl = jpl 451 451 DO_2D( 1, 1, 1, 1 ) 452 452 IF( h_i(ji,jj,jl) > 50._wp ) THEN … … 460 460 ialert_id = ialert_id + 1 ! reference number of this alert 461 461 cl_alname(ialert_id) = ' Very thin ice ' ! name of the alert 462 jl = 1 462 jl = 1 463 463 DO_2D( 1, 1, 1, 1 ) 464 464 IF( h_i(ji,jj,jl) < rn_himin ) THEN … … 484 484 cl_alname(ialert_id) = ' Ice on continents ' ! name of the alert 485 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 486 IF( tmask(ji,jj,1) == 0._wp .AND. ( at_i(ji,jj) > 0._wp .OR. vt_i(ji,jj) > 0._wp ) ) THEN 487 487 WRITE(numout,*) ' ALERTE : Ice on continents ',at_i(ji,jj),vt_i(ji,jj) 488 488 WRITE(numout,*) ' at i,j = ',ji,jj … … 496 496 DO_2D( 1, 1, 1, 1 ) 497 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 498 & ( vt_i(ji,jj) > 0._wp .AND. at_i(ji,jj) == 0._wp ) ) THEN 499 499 WRITE(numout,*) ' ALERTE : Incompatible ice conc and vol ',at_i(ji,jj),vt_i(ji,jj) 500 500 WRITE(numout,*) ' at i,j = ',ji,jj … … 520 520 ! 521 521 END SUBROUTINE ice_ctl 522 522 523 523 SUBROUTINE ice_prt( kt, ki, kj, kn, cd1 ) 524 524 !!------------------------------------------------------------------- 525 !! *** ROUTINE ice_prt *** 526 !! 527 !! ** Purpose : Writes global ice state on the (i,j) point 528 !! in ocean.ouput 529 !! 3 possibilities exist 525 !! *** ROUTINE ice_prt *** 526 !! 527 !! ** Purpose : Writes global ice state on the (i,j) point 528 !! in ocean.ouput 529 !! 3 possibilities exist 530 530 !! n = 1/-1 -> simple ice state 531 531 !! n = 2 -> exhaustive state 532 532 !! n = 3 -> ice/ocean salt fluxes 533 533 !! 534 !! ** input : point coordinates (i,j) 534 !! ** input : point coordinates (i,j) 535 535 !! n : number of the option 536 536 !!------------------------------------------------------------------- … … 550 550 ! Simple state 551 551 !---------------- 552 552 553 553 IF ( kn == 1 .OR. kn == -1 ) THEN 554 554 WRITE(numout,*) ' ice_prt - Point : ',ji,jj … … 566 566 WRITE(numout,*) ' - Cell values ' 567 567 WRITE(numout,*) ' ~~~~~~~~~~~ ' 568 WRITE(numout,*) ' at_i : ', at_i(ji,jj) 569 WRITE(numout,*) ' ato_i : ', ato_i(ji,jj) 570 WRITE(numout,*) ' vt_i : ', vt_i(ji,jj) 571 WRITE(numout,*) ' vt_s : ', vt_s(ji,jj) 568 WRITE(numout,*) ' at_i : ', at_i(ji,jj) 569 WRITE(numout,*) ' ato_i : ', ato_i(ji,jj) 570 WRITE(numout,*) ' vt_i : ', vt_i(ji,jj) 571 WRITE(numout,*) ' vt_s : ', vt_s(ji,jj) 572 572 DO jl = 1, jpl 573 573 WRITE(numout,*) ' - Category (', jl,')' … … 592 592 ! Exhaustive state 593 593 !-------------------- 594 594 595 595 IF ( kn .EQ. 2 ) THEN 596 596 WRITE(numout,*) ' ice_prt - Point : ',ji,jj … … 598 598 WRITE(numout,*) ' Exhaustive state ' 599 599 WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 600 WRITE(numout,*) 600 WRITE(numout,*) 601 601 WRITE(numout,*) ' - Cell values ' 602 602 WRITE(numout,*) ' ~~~~~~~~~~~ ' 603 WRITE(numout,*) ' at_i : ', at_i(ji,jj) 604 WRITE(numout,*) ' vt_i : ', vt_i(ji,jj) 605 WRITE(numout,*) ' vt_s : ', vt_s(ji,jj) 603 WRITE(numout,*) ' at_i : ', at_i(ji,jj) 604 WRITE(numout,*) ' vt_i : ', vt_i(ji,jj) 605 WRITE(numout,*) ' vt_s : ', vt_s(ji,jj) 606 606 WRITE(numout,*) ' u_ice(i-1,j) : ', u_ice(ji-1,jj) 607 607 WRITE(numout,*) ' u_ice(i ,j) : ', u_ice(ji,jj) … … 610 610 WRITE(numout,*) ' strength : ', strength(ji,jj) 611 611 WRITE(numout,*) 612 612 613 613 DO jl = 1, jpl 614 614 WRITE(numout,*) ' - Category (',jl,')' 615 WRITE(numout,*) ' ~~~~~~~~ ' 615 WRITE(numout,*) ' ~~~~~~~~ ' 616 616 WRITE(numout,*) ' h_i : ', h_i(ji,jj,jl) , ' h_s : ', h_s(ji,jj,jl) 617 617 WRITE(numout,*) ' t_i : ', t_i(ji,jj,1:nlay_i,jl) 618 618 WRITE(numout,*) ' t_su : ', t_su(ji,jj,jl) , ' t_s : ', t_s(ji,jj,1:nlay_s,jl) 619 619 WRITE(numout,*) ' s_i : ', s_i(ji,jj,jl) , ' o_i : ', o_i(ji,jj,jl) 620 WRITE(numout,*) ' a_i : ', a_i(ji,jj,jl) , ' a_i_b : ', a_i_b(ji,jj,jl) 621 WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl) , ' v_i_b : ', v_i_b(ji,jj,jl) 622 WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl) , ' v_s_b : ', v_s_b(ji,jj,jl) 623 WRITE(numout,*) ' e_i1 : ', e_i(ji,jj,1,jl) , ' ei1 : ', e_i_b(ji,jj,1,jl) 624 WRITE(numout,*) ' e_i2 : ', e_i(ji,jj,2,jl) , ' ei2_b : ', e_i_b(ji,jj,2,jl) 625 WRITE(numout,*) ' e_snow : ', e_s(ji,jj,1,jl) , ' e_snow_b : ', e_s_b(ji,jj,1,jl) 626 WRITE(numout,*) ' sv_i : ', sv_i(ji,jj,jl) , ' sv_i_b : ', sv_i_b(ji,jj,jl) 620 WRITE(numout,*) ' a_i : ', a_i(ji,jj,jl) , ' a_i_b : ', a_i_b(ji,jj,jl) 621 WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl) , ' v_i_b : ', v_i_b(ji,jj,jl) 622 WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl) , ' v_s_b : ', v_s_b(ji,jj,jl) 623 WRITE(numout,*) ' e_i1 : ', e_i(ji,jj,1,jl) , ' ei1 : ', e_i_b(ji,jj,1,jl) 624 WRITE(numout,*) ' e_i2 : ', e_i(ji,jj,2,jl) , ' ei2_b : ', e_i_b(ji,jj,2,jl) 625 WRITE(numout,*) ' e_snow : ', e_s(ji,jj,1,jl) , ' e_snow_b : ', e_s_b(ji,jj,1,jl) 626 WRITE(numout,*) ' sv_i : ', sv_i(ji,jj,jl) , ' sv_i_b : ', sv_i_b(ji,jj,jl) 627 627 END DO !jl 628 628 629 629 WRITE(numout,*) 630 630 WRITE(numout,*) ' - Heat / FW fluxes ' … … 634 634 WRITE(numout,*) ' qns_ini : ', (1._wp-at_i_b(ji,jj)) * qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) ) 635 635 WRITE(numout,*) 636 WRITE(numout,*) 637 WRITE(numout,*) ' sst : ', sst_m(ji,jj) 638 WRITE(numout,*) ' sss : ', sss_m(ji,jj) 639 WRITE(numout,*) 636 WRITE(numout,*) 637 WRITE(numout,*) ' sst : ', sst_m(ji,jj) 638 WRITE(numout,*) ' sss : ', sss_m(ji,jj) 639 WRITE(numout,*) 640 640 WRITE(numout,*) ' - Stresses ' 641 641 WRITE(numout,*) ' ~~~~~~~~ ' 642 WRITE(numout,*) ' utau_ice : ', utau_ice(ji,jj) 642 WRITE(numout,*) ' utau_ice : ', utau_ice(ji,jj) 643 643 WRITE(numout,*) ' vtau_ice : ', vtau_ice(ji,jj) 644 WRITE(numout,*) ' utau : ', utau (ji,jj) 644 WRITE(numout,*) ' utau : ', utau (ji,jj) 645 645 WRITE(numout,*) ' vtau : ', vtau (ji,jj) 646 646 ENDIF 647 647 648 648 !--------------------- 649 649 ! Salt / heat fluxes 650 650 !--------------------- 651 651 652 652 IF ( kn .EQ. 3 ) THEN 653 653 WRITE(numout,*) ' ice_prt - Point : ',ji,jj … … 664 664 WRITE(numout,*) ' qt_atm_oi : ', qt_atm_oi(ji,jj) 665 665 WRITE(numout,*) ' qt_oce_ai : ', qt_oce_ai(ji,jj) 666 WRITE(numout,*) ' dhc : ', diag_heat(ji,jj) 666 WRITE(numout,*) ' dhc : ', diag_heat(ji,jj) 667 667 WRITE(numout,*) 668 668 WRITE(numout,*) ' hfx_dyn : ', hfx_dyn(ji,jj) 669 669 WRITE(numout,*) ' hfx_thd : ', hfx_thd(ji,jj) 670 670 WRITE(numout,*) ' hfx_res : ', hfx_res(ji,jj) 671 WRITE(numout,*) ' qsb_ice_bot : ', qsb_ice_bot(ji,jj) 671 WRITE(numout,*) ' qsb_ice_bot : ', qsb_ice_bot(ji,jj) 672 672 WRITE(numout,*) ' qlead : ', qlead(ji,jj) * r1_Dt_ice 673 673 WRITE(numout,*) … … 680 680 WRITE(numout,*) 681 681 WRITE(numout,*) ' - Momentum fluxes ' 682 WRITE(numout,*) ' utau : ', utau(ji,jj) 682 WRITE(numout,*) ' utau : ', utau(ji,jj) 683 683 WRITE(numout,*) ' vtau : ', vtau(ji,jj) 684 ENDIF 684 ENDIF 685 685 WRITE(numout,*) ' ' 686 686 ! … … 694 694 !! *** ROUTINE ice_prt3D *** 695 695 !! 696 !! ** Purpose : CTL prints of ice arrays in case sn_cfctl%prtctl is activated 696 !! ** Purpose : CTL prints of ice arrays in case sn_cfctl%prtctl is activated 697 697 !! 698 698 !!------------------------------------------------------------------- 699 699 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 700 700 INTEGER :: jk, jl ! dummy loop indices 701 701 702 702 CALL prt_ctl_info(' ========== ') 703 703 CALL prt_ctl_info( cd_routine ) … … 718 718 CALL prt_ctl(tab2d_1=delta_i , clinfo1=' delta_i :') 719 719 CALL prt_ctl(tab2d_1=u_ice , clinfo1=' u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 720 720 721 721 DO jl = 1, jpl 722 722 CALL prt_ctl_info(' ') … … 735 735 CALL prt_ctl(tab2d_1=sv_i (:,:,jl) , clinfo1= ' sv_i : ') 736 736 CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' oa_i : ') 737 737 738 738 DO jk = 1, nlay_i 739 739 CALL prt_ctl_info(' - Layer : ', ivar=jk) … … 742 742 END DO 743 743 END DO 744 744 745 745 CALL prt_ctl_info(' ') 746 746 CALL prt_ctl_info(' - Stresses : ') … … 748 748 CALL prt_ctl(tab2d_1=utau , clinfo1= ' utau : ', tab2d_2=vtau , clinfo2= ' vtau : ') 749 749 CALL prt_ctl(tab2d_1=utau_ice , clinfo1= ' utau_ice : ', tab2d_2=vtau_ice , clinfo2= ' vtau_ice : ') 750 750 751 751 END SUBROUTINE ice_prt3D 752 752 … … 853 853 !!---------------------------------------------------------------------- 854 854 !! *** ROUTINE ice_drift_init *** 855 !! 855 !! 856 856 !! ** Purpose : create output file, initialise arrays 857 857 !!---------------------------------------------------------------------- … … 879 879 ! 880 880 END SUBROUTINE ice_drift_init 881 881 882 882 #else 883 883 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
