Changeset 13601 for NEMO/trunk
- Timestamp:
- 2020-10-14T17:59:34+02:00 (4 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 11 edited
-
cfgs/SHARED/field_def_nemo-ice.xml (modified) (3 diffs)
-
src/ICE/ice.F90 (modified) (2 diffs)
-
src/ICE/icecor.F90 (modified) (1 diff)
-
src/ICE/icectl.F90 (modified) (7 diffs)
-
src/ICE/icedyn_rhg_evp.F90 (modified) (5 diffs)
-
src/ICE/icestp.F90 (modified) (4 diffs)
-
src/ICE/icethd.F90 (modified) (10 diffs)
-
src/ICE/icethd_dh.F90 (modified) (1 diff)
-
src/ICE/icethd_do.F90 (modified) (3 diffs)
-
src/ICE/iceupdate.F90 (modified) (8 diffs)
-
src/OCE/BDY/bdyice.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/cfgs/SHARED/field_def_nemo-ice.xml
r13472 r13601 78 78 <field id="isig2" long_name="2nd principal stress component for EVP rhg" unit="" /> 79 79 <field id="isig3" long_name="convergence measure for EVP rheology (must be around 1)" unit="" /> 80 <!-- clem: sig1 sig2 sig 3 will be replaced by sig1 and sig2 in a following commit --> 81 <field id="sig1_pnorm" long_name="P-normalized 1st principal stress component" unit="" /> 82 <field id="sig2_pnorm" long_name="P-normalized 2nd principal stress component" unit="" /> 80 83 <field id="normstr" long_name="Average normal stress in sea ice" standard_name="average_normal_stress" unit="N/m" /> 81 84 <field id="sheastr" long_name="Maximum shear stress in sea ice" standard_name="maximum_shear_stress" unit="N/m" /> … … 165 168 166 169 <!-- diags --> 167 <field id="hfxdhc" long_name="Heat content variation in snow and ice (neg = ice cooling)" unit="W/m2" /> 170 <field id="hfxdhc" long_name="Heat content variation in snow and ice (neg = ice cooling)" unit="W/m2" /> 171 <!-- available if ln_icediachk=T --> 172 <field id="icedrift_mass" long_name="Ice mass drift (conservation check)" unit="kg/m2/s" /> 173 <field id="icedrift_salt" long_name="Ice salt drift (conservation check)" unit="kg/m2/s" /> 174 <field id="icedrift_heat" long_name="Ice heat drift (conservation check)" unit="W/m2" /> 168 175 169 176 <!-- sbcssm variables --> … … 459 466 <field field_ref="vfxthin" name="vfxthin" /> 460 467 461 <!-- diag error for negative ice volume after advection -->462 <field field_ref="iceneg_pres" name="sineg_pres" />463 <field field_ref="iceneg_volu" name="sineg_volu" />464 <field field_ref="iceneg_hfx" name="sineg_hfx" />465 468 </field_group> 466 469 -
NEMO/trunk/src/ICE/ice.F90
r13472 r13601 392 392 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 393 393 ! 394 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_mass !: advection of mass (kg/m2/s) 395 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_salt !: advection of salt (g/m2/s) 396 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_heat !: advection of heat (W/m2) 397 ! 394 398 !!---------------------------------------------------------------------- 395 399 !! * Ice conservation … … 495 499 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & 496 500 & diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), & 497 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), STAT=ierr(ii) ) 501 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), & 502 & diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) 498 503 499 504 ! * Ice conservation -
NEMO/trunk/src/ICE/icecor.F90
r13497 r13601 95 95 zsal = sv_i(ji,jj,jl) 96 96 sv_i(ji,jj,jl) = MIN( MAX( rn_simin*v_i(ji,jj,jl) , sv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl) ) 97 sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc ! associated salt flux 97 IF( kn /= 0 ) & ! no ice-ocean exchanges if kn=0 (for bdy for instance) otherwise conservation diags will fail 98 & sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc ! associated salt flux 98 99 END_2D 99 100 END DO 100 101 ENDIF 101 ! !-----------------------------------------------------102 CALL ice_var_zapsmall ! Zap small values !103 ! !-----------------------------------------------------104 102 103 IF( kn /= 0 ) THEN ! no zapsmall if kn=0 (for bdy for instance) because we do not want ice-ocean exchanges (wfx,sfx,hfx) 104 ! otherwise conservation diags will fail 105 ! !----------------------------------------------------- 106 CALL ice_var_zapsmall ! Zap small values ! 107 ! !----------------------------------------------------- 108 ENDIF 105 109 ! !----------------------------------------------------- 106 110 IF( kn == 2 ) THEN ! Ice drift case: Corrections to avoid wrong values ! -
NEMO/trunk/src/ICE/icectl.F90
r13472 r13601 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 … … 132 140 133 141 ! -- 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)142 zvtrp = glob_sum( 'icectl', diag_adv_mass * e1e2t ) 143 zetrp = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 136 144 137 145 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 156 164 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_amax 157 165 ! 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 166 IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 167 & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 168 IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 169 & WRITE(numout,*) cd_routine,' : violation adv scheme [J] = ',zetrp * rdt_ice 163 170 ENDIF 164 171 ! … … 186 193 ! water flux 187 194 ! -- mass diag -- ! 188 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) 195 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub & 196 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) 189 197 190 198 ! -- salt diag -- ! 191 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t )199 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) 192 200 193 201 ! -- 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 ) 202 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 203 ! equivalent to this: 204 !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 205 !! & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr & 206 !! & ) * e1e2t ) 197 207 198 208 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 204 214 IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 205 215 & 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 216 IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 217 & WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 207 218 ENDIF 208 219 ! … … 725 736 726 737 END SUBROUTINE ice_prt3D 738 739 740 SUBROUTINE ice_drift_wri( kt ) 741 !!------------------------------------------------------------------- 742 !! *** ROUTINE ice_drift_wri *** 743 !! 744 !! ** Purpose : conservation of mass, salt and heat 745 !! write the drift in a ascii file at each time step 746 !! and the total run drifts 747 !!------------------------------------------------------------------- 748 INTEGER, INTENT(in) :: kt ! ice time-step index 749 ! 750 INTEGER :: ji, jj 751 REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 752 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_mass2D, zdiag_salt2D, zdiag_heat2D 753 !!------------------------------------------------------------------- 754 ! 755 IF( kt == nit000 .AND. lwp ) THEN 756 WRITE(numout,*) 757 WRITE(numout,*) 'ice_drift_wri: sea-ice drifts' 758 WRITE(numout,*) '~~~~~~~~~~~~~' 759 ENDIF 760 ! 761 ! 2D budgets (must be close to 0) 762 IF( iom_use('icedrift_mass') .OR. iom_use('icedrift_salt') .OR. iom_use('icedrift_heat') ) THEN 763 DO_2D( 1, 1, 1, 1 ) 764 zdiag_mass2D(ji,jj) = wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_spr(ji,jj) + wfx_sub(ji,jj) & 765 & + diag_vice(ji,jj) + diag_vsnw(ji,jj) - diag_adv_mass(ji,jj) 766 zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) 767 zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj) 768 END_2D 769 ! 770 ! write outputs 771 CALL iom_put( 'icedrift_mass', zdiag_mass2D ) 772 CALL iom_put( 'icedrift_salt', zdiag_salt2D ) 773 CALL iom_put( 'icedrift_heat', zdiag_heat2D ) 774 ENDIF 775 776 ! -- mass diag -- ! 777 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub & 778 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rdt_ice 779 zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 780 781 ! -- salt diag -- ! 782 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rdt_ice * 1.e-3 783 zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 784 785 ! -- heat diag -- ! 786 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 787 zdiag_adv_heat = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 788 789 ! ! write out to file 790 IF( lwp ) THEN 791 ! check global drift (must be close to 0) 792 WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift [kg]', zdiag_mass 793 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift [kg]', zdiag_salt 794 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift [W] ', zdiag_heat 795 ! check drift from advection scheme (can be /=0 with bdy but not sure why) 796 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'mass drift adv [kg]', zdiag_adv_mass 797 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift adv [kg]', zdiag_adv_salt 798 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift adv [W] ', zdiag_adv_heat 799 ENDIF 800 ! ! drifts 801 rdiag_icemass = rdiag_icemass + zdiag_mass 802 rdiag_icesalt = rdiag_icesalt + zdiag_salt 803 rdiag_iceheat = rdiag_iceheat + zdiag_heat 804 rdiag_adv_icemass = rdiag_adv_icemass + zdiag_adv_mass 805 rdiag_adv_icesalt = rdiag_adv_icesalt + zdiag_adv_salt 806 rdiag_adv_iceheat = rdiag_adv_iceheat + zdiag_adv_heat 807 ! 808 ! ! output drifts and close ascii file 809 IF( kt == nitend - nn_fsbc + 1 .AND. lwp ) THEN 810 ! to ascii file 811 WRITE(numicedrift,*) '******************************************' 812 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift [kg]', rdiag_icemass 813 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift adv [kg]', rdiag_adv_icemass 814 WRITE(numicedrift,*) '******************************************' 815 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift [kg]', rdiag_icesalt 816 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift adv [kg]', rdiag_adv_icesalt 817 WRITE(numicedrift,*) '******************************************' 818 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift [W] ', rdiag_iceheat 819 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift adv [W] ', rdiag_adv_iceheat 820 CLOSE( numicedrift ) 821 ! 822 ! to ocean output 823 WRITE(numout,*) 824 WRITE(numout,*) 'ice_drift_wri: ice drifts information for the run ' 825 WRITE(numout,*) '~~~~~~~~~~~~~' 826 ! check global drift (must be close to 0) 827 WRITE(numout,*) ' sea-ice mass drift [kg] = ', rdiag_icemass 828 WRITE(numout,*) ' sea-ice salt drift [kg] = ', rdiag_icesalt 829 WRITE(numout,*) ' sea-ice heat drift [W] = ', rdiag_iceheat 830 ! check drift from advection scheme (can be /=0 with bdy but not sure why) 831 WRITE(numout,*) ' sea-ice mass drift adv [kg] = ', rdiag_adv_icemass 832 WRITE(numout,*) ' sea-ice salt drift adv [kg] = ', rdiag_adv_icesalt 833 WRITE(numout,*) ' sea-ice heat drift adv [W] = ', rdiag_adv_iceheat 834 ENDIF 835 ! 836 END SUBROUTINE ice_drift_wri 837 838 SUBROUTINE ice_drift_init 839 !!---------------------------------------------------------------------- 840 !! *** ROUTINE ice_drift_init *** 841 !! 842 !! ** Purpose : create output file, initialise arrays 843 !!---------------------------------------------------------------------- 844 ! 845 IF( .NOT.ln_icediachk ) RETURN ! exit 846 ! 847 IF(lwp) THEN 848 WRITE(numout,*) 849 WRITE(numout,*) 'ice_drift_init: Output ice drifts to ',TRIM(clname), ' file' 850 WRITE(numout,*) '~~~~~~~~~~~~~' 851 WRITE(numout,*) 852 ! 853 ! create output ascii file 854 CALL ctl_opn( numicedrift, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, narea ) 855 WRITE(numicedrift,*) 'Timestep Drifts' 856 WRITE(numicedrift,*) '******************************************' 857 ENDIF 858 ! 859 rdiag_icemass = 0._wp 860 rdiag_icesalt = 0._wp 861 rdiag_iceheat = 0._wp 862 rdiag_adv_icemass = 0._wp 863 rdiag_adv_icesalt = 0._wp 864 rdiag_adv_iceheat = 0._wp 865 ! 866 END SUBROUTINE ice_drift_init 727 867 728 868 #else -
NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
r13550 r13601 170 170 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 171 171 !! --- diags 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 172 REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength 173 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p, zten_i 173 174 !! --- SIMIP diags 174 175 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 726 727 & ) * r1_e1e2t(ji,jj) 727 728 zdt2 = zdt * zdt 729 730 zten_i(ji,jj) = zdt 728 731 729 732 ! shear**2 at T points (doc eq. A16) … … 741 744 742 745 ! delta at T points 743 z delta(ji,jj) = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )744 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -z delta(ji,jj)) ) ! 0 if delta=0745 pdelta_i(ji,jj) = z delta(ji,jj) + rn_creepl * rswitch746 zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 747 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 748 pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 746 749 747 750 END_2D 748 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp ) 751 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, & 752 & zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp ) 749 753 750 754 ! --- Store the stress tensor for the next time step --- ! 751 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )752 755 pstress1_i (:,:) = zs1 (:,:) 753 756 pstress2_i (:,:) = zs2 (:,:) … … 778 781 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 779 782 780 ! --- stress tensor--- !781 IF( iom_use(' isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN782 ! 783 ALLOCATE( zsig 1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )783 ! --- Stress tensor invariants (SIMIP diags) --- ! 784 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 785 ! 786 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 784 787 ! 785 DO_2D( 0, 0, 0, 0 ) 786 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 787 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 788 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 789 790 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 791 792 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 793 794 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 795 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 796 !! zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 797 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 798 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 799 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 800 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 801 END_2D 802 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp ) 803 ! 804 CALL iom_put( 'isig1' , zsig1 ) 805 CALL iom_put( 'isig2' , zsig2 ) 806 CALL iom_put( 'isig3' , zsig3 ) 807 ! 808 ! Stress tensor invariants (normal and shear stress N/m) 809 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 810 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 811 812 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 788 DO_2D( 1, 1, 1, 1 ) 789 790 ! Ice stresses 791 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 792 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 793 ! I know, this can be confusing... 794 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 795 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 796 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 797 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 798 799 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 800 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 801 zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 ) ! 2nd '' '', aka maximum shear stress 802 803 END_2D 804 ! 805 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 806 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 807 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 808 809 DEALLOCATE ( zsig_I, zsig_II ) 810 811 ENDIF 812 813 ! --- Normalized stress tensor principal components --- ! 814 ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 815 ! Recommendation 1 : we use ice strength, not replacement pressure 816 ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 817 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 818 ! 819 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 820 ! 821 DO_2D( 1, 1, 1, 1 ) 822 823 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 824 ! and **deformations** at current iterates 825 ! following Lemieux & Dupont (2020) 826 zfac = zp_delt(ji,jj) 827 zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 828 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 829 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 830 831 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 832 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 833 zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 ) ! 2nd '' '', aka maximum shear stress 834 835 ! Normalized principal stresses (used to display the ellipse) 836 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 837 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 838 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 839 END_2D 840 ! 841 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 842 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 843 844 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 845 813 846 ENDIF 814 847 … … 946 979 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 947 980 ! close file 948 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid)981 IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid) 949 982 ENDIF 950 983 -
NEMO/trunk/src/ICE/icestp.F90
r13558 r13601 197 197 IF( ln_icediahsb ) CALL ice_dia( kt ) ! -- Diagnostics outputs 198 198 ! 199 IF( ln_icediachk ) CALL ice_drift_wri( kt ) ! -- Diagnostics outputs for conservation 200 ! 199 201 CALL ice_wri( kt ) ! -- Ice outputs 200 202 ! … … 281 283 ! 282 284 CALL ice_dia_init ! initialization for diags 285 ! 286 CALL ice_drift_init ! initialization for diags of conservation 283 287 ! 284 288 fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction … … 341 345 ENDIF 342 346 ! 343 IF( ln_bdy .AND. ln_icediachk ) CALL ctl_warn('par_init: online conservation check does not work with BDY')344 !345 347 rDt_ice = REAL(nn_fsbc) * rn_Dt !--- sea-ice timestep and its inverse 346 348 r1_Dt_ice = 1._wp / rDt_ice … … 438 440 diag_trp_ei(:,:) = 0._wp ; diag_trp_es(:,:) = 0._wp 439 441 diag_trp_sv(:,:) = 0._wp 442 diag_adv_mass(:,:) = 0._wp 443 diag_adv_salt(:,:) = 0._wp 444 diag_adv_heat(:,:) = 0._wp 440 445 441 446 END SUBROUTINE diag_set0 -
NEMO/trunk/src/ICE/icethd.F90
r13547 r13601 18 18 USE ice ! sea-ice: variables 19 19 !!gm list trop longue ==>>> why not passage en argument d'appel ? 20 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot,sprecip, ln_cpl20 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl 21 21 USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 22 22 & qml_ice, qcn_ice, qtr_ice_top … … 52 52 LOGICAL :: ln_icedO ! activate ice growth in open-water (T) or not (F) 53 53 LOGICAL :: ln_icedS ! activate gravity drainage and flushing (T) or not (F) 54 LOGICAL :: ln_leadhfx ! 54 LOGICAL :: ln_leadhfx ! heat in the leads is used to melt sea-ice before warming the ocean 55 55 56 56 !! for convergence tests … … 91 91 ! 92 92 INTEGER :: ji, jj, jk, jl ! dummy loop indices 93 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg 94 REAL(wp), PARAMETER :: zfric_umin = 0._wp 95 REAL(wp), PARAMETER :: zch = 0.0057_wp 96 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)93 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos 94 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 95 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 96 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 97 97 ! 98 98 !!------------------------------------------------------------------- … … 124 124 & ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 125 125 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1) 126 zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + & 127 & ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 126 128 END_2D 127 129 ELSE ! if no ice dynamics => transmit directly the atmospheric stress to the ocean … … 130 132 & ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 131 133 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 134 zvel(ji,jj) = 0._wp 132 135 END_2D 133 136 ENDIF 134 CALL lbc_lnk ( 'icethd', zfric, 'T',1.0_wp )137 CALL lbc_lnk_multi( 'icethd', zfric, 'T', 1.0_wp, zvel, 'T', 1.0_wp ) 135 138 ! 136 139 !--------------------------------------------------------------------! … … 140 143 rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 141 144 ! 142 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget143 ! ! practically no "direct lateral ablation"144 !145 ! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean146 ! ! temperature and turbulent mixing (McPhee, 1992)147 !148 145 ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 149 146 zqld = tmask(ji,jj,1) * rDt_ice * & … … 151 148 & ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 152 149 153 ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- ! 150 ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 151 ! (mostly<0 but >0 if supercooling) 154 152 zqfr = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1) ! both < 0 (t_bo < sst) and > 0 (t_bo > sst) 155 153 zqfr_neg = MIN( zqfr , 0._wp ) ! only < 0 156 157 ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2) 154 zqfr_pos = MAX( zqfr , 0._wp ) ! only > 0 155 156 ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 157 ! (mostly>0 but <0 if supercooling) 158 158 zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 159 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 160 161 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 159 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 160 162 161 ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach 163 162 ! the freezing point, so that we do not have SST < T_freeze 164 ! This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0 165 166 !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice 167 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 168 169 ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting 170 ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 171 IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 172 IF( ln_leadhfx ) THEN ; fhld(ji,jj) = rswitch * zqld * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 173 ELSE ; fhld(ji,jj) = 0._wp 174 ENDIF 163 ! This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 164 ! The following formulation is ok for both normal conditions and supercooling 165 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 166 167 ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 168 ! qlead is the energy received from the atm. in the leads. 169 ! If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld (W/m2) 170 ! If cooling (zqld < 0), then the energy in the leads is used to grow ice in open water => qlead (J.m-2) 171 IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 172 ! upper bound for fhld: fhld should be equal to zqld 173 ! but we have to make sure that this heat will not make the sst drop below the freezing point 174 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 175 ! The following formulation is ok for both normal conditions and supercooling 176 fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) & ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 177 & - qsb_ice_bot(ji,jj) ) 175 178 qlead(ji,jj) = 0._wp 176 179 ELSE 177 180 fhld (ji,jj) = 0._wp 181 ! upper bound for qlead: qlead should be equal to zqld 182 ! but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 183 ! The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 184 ! and freezing point is reached if zqfr = zqld - qsb*a/dt 185 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 186 ! The following formulation is ok for both normal conditions and supercooling 187 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 178 188 ENDIF 179 189 ! 180 ! Net heat flux on top of the ice-ocean [W.m-2] 181 ! --------------------------------------------- 182 qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 190 ! If ice is landfast and ice concentration reaches its max 191 ! => stop ice formation in open water 192 IF( zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 ) qlead(ji,jj) = 0._wp 193 ! 194 ! If the grid cell is almost fully covered by ice (no leads) 195 ! => stop ice formation in open water 196 IF( at_i(ji,jj) >= (1._wp - epsi10) ) qlead(ji,jj) = 0._wp 197 ! 198 ! If ln_leadhfx is false 199 ! => do not use energy of the leads to melt sea-ice 200 IF( .NOT.ln_leadhfx ) fhld(ji,jj) = 0._wp 201 ! 183 202 END_2D 184 203 … … 191 210 ENDIF 192 211 193 ! ---------------------------------------------------------------------194 ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]195 ! ---------------------------------------------------------------------196 ! First step here : non solar + precip - qlead - qsensible197 ! Second step in icethd_dh : heat remaining if total melt (zq_rema)198 ! Third step in iceupdate.F90 : heat from ice-ocean mass exchange (zf_mass) + solar199 qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) & ! Non solar heat flux received by the ocean200 & - qlead(:,:) * r1_Dt_ice & ! heat flux taken from the ocean where there is open water ice formation201 & - at_i (:,:) * qsb_ice_bot(:,:) & ! heat flux taken by sensible flux202 & - at_i (:,:) * fhld (:,:) ! heat flux taken during bottom growth/melt203 ! ! (fhld should be 0 while bott growth)204 212 !-------------------------------------------------------------------------------------------! 205 213 ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories … … 418 426 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 419 427 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 420 CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )421 428 ! 422 429 ! ocean surface fields … … 510 517 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 511 518 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 512 CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )513 519 ! 514 520 CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d (1:npti), qns_ice (:,:,kl) ) -
NEMO/trunk/src/ICE/icethd_dh.F90
r13472 r13601 556 556 ! 557 557 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 558 qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice558 !!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 559 559 560 560 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) -
NEMO/trunk/src/ICE/icethd_do.F90
r13547 r13601 131 131 132 132 ! Default new ice thickness 133 WHERE( qlead(:,:) < 0._wp .AND. tau_icebfr(:,:) == 0._wp ) ; ht_i_new(:,:) = rn_hinew ! if cooling and no landfast 134 ELSEWHERE ; ht_i_new(:,:) = 0._wp 133 WHERE( qlead(:,:) < 0._wp ) ! cooling 134 ht_i_new(:,:) = rn_hinew 135 ELSEWHERE 136 ht_i_new(:,:) = 0._wp 135 137 END WHERE 136 138 … … 146 148 ! 147 149 DO_2D( 0, 0, 0, 0 ) 148 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast150 IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 149 151 ! -- Wind stress -- ! 150 152 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & … … 198 200 ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice 199 201 !------------------------------------------------------------------------------! 200 ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice202 ! it occurs if cooling 201 203 202 204 ! Identify grid points where new ice forms 203 205 npti = 0 ; nptidx(:) = 0 204 206 DO_2D( 1, 1, 1, 1 ) 205 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp) THEN207 IF ( qlead(ji,jj) < 0._wp ) THEN 206 208 npti = npti + 1 207 209 nptidx( npti ) = (jj - 1) * jpi + ji -
NEMO/trunk/src/ICE/iceupdate.F90
r13497 r13601 24 24 USE traqsr ! add penetration of solar flux in the calculation of heat budget 25 25 USE icectl ! sea-ice: control prints 26 USE bdy_oce , ONLY : ln_bdy27 26 USE zdfdrg , ONLY : ln_drgice_imp 28 27 ! … … 92 91 ! 93 92 INTEGER :: ji, jj, jl, jk ! dummy loop indices 94 REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2)95 93 REAL(wp) :: zqsr ! New solar flux received by the ocean 96 94 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace … … 103 101 WRITE(numout,*)'~~~~~~~~~~~~~~' 104 102 ENDIF 103 104 ! Net heat flux on top of the ice-ocean (W.m-2) 105 !---------------------------------------------- 106 qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:) 105 107 106 108 ! --- case we bypass ice thermodynamics --- ! … … 115 117 DO_2D( 1, 1, 1, 1 ) 116 118 117 ! Solar heat flux reaching the ocean = zqsr (W.m-2)119 ! Solar heat flux reaching the ocean (max) = zqsr (W.m-2) 118 120 !--------------------------------------------------- 119 121 zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) ) … … 121 123 ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2) 122 124 !--------------------------------------------------- 123 zqmass = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 124 qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + zqmass + zqsr 125 126 ! Add the residual from heat diffusion equation and sublimation (W.m-2) 127 !---------------------------------------------------------------------- 128 qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + hfx_err_dif(ji,jj) + & 129 & ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) ) 130 125 qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) & 126 & - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) & 127 & + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) & 128 & + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj) 129 131 130 ! New qsr and qns used to compute the oceanic heat flux at the next time step 132 131 !---------------------------------------------------------------------------- 133 qsr(ji,jj) = zqsr 132 ! if warming and some ice remains, then we suppose that the whole solar flux has been consumed to melt the ice 133 ! else ( cooling or no ice left ), then we suppose that no solar flux has been consumed 134 ! 135 IF( fhld(ji,jj) > 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN !-- warming and some ice remains 136 ! solar flux transmitted thru the 1st level of the ocean (i.e. not used by sea-ice) 137 qsr(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * ( 1._wp - frq_m(ji,jj) ) & 138 ! + solar flux transmitted thru ice (also not used by sea-ice) 139 & + SUM( a_i_b(ji,jj,:) * qtr_ice_bot(ji,jj,:) ) 140 ! 141 ELSE !-- cooling or no ice left 142 qsr(ji,jj) = zqsr 143 ENDIF 144 ! 145 ! the non-solar is simply derived from the solar flux 134 146 qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr 135 147 136 148 ! Mass flux at the atm. surface 137 149 !----------------------------------- … … 140 152 ! Mass flux at the ocean surface 141 153 !------------------------------------ 142 ! case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED) 143 ! ------------------------------------------------------------------------------------- 144 ! The idea of this approach is that the system that we consider is the ICE-OCEAN system 145 ! Thus FW flux = External ( E-P+snow melt) 146 ! Salt flux = Exchanges in the ice-ocean system then converted into FW 147 ! Associated to Ice formation AND Ice melting 148 ! Even if i see Ice melting as a FW and SALT flux 149 ! 150 ! mass flux from ice/ocean 154 ! ice-ocean mass flux 151 155 wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) & 152 156 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 153 154 ! add the snow melt water to snow mass flux to the ocean157 158 ! snw-ocean mass flux 155 159 wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj) 156 157 ! mass flux at the ocean/ice interface 158 fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) ) ! F/M mass flux save at least for biogeochemical model 159 emp(ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! mass flux + F/M mass flux (always ice/ocean mass exchange) 160 160 161 ! total mass flux at the ocean/ice interface 162 fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model 163 emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux 161 164 162 165 ! Salt flux at the ocean surface … … 262 265 CALL iom_put ('hfxdif' , hfx_dif ) ! heat flux used for ice temperature change 263 266 CALL iom_put ('hfxsnw' , hfx_snw ) ! heat flux used for snow melt 264 CALL iom_put ('hfxerr' , hfx_err_dif ) ! heat flux error after heat diffusion (included in qt_oce_ai)267 CALL iom_put ('hfxerr' , hfx_err_dif ) ! heat flux error after heat diffusion 265 268 266 269 ! heat fluxes associated with mass exchange (freeze/melt/precip...) … … 279 282 !--------- 280 283 #if ! defined key_agrif 281 IF( ln_icediachk .AND. .NOT. ln_bdy) CALL ice_cons_final('iceupdate') ! conservation284 IF( ln_icediachk ) CALL ice_cons_final('iceupdate') ! conservation 282 285 #endif 283 IF( ln_icectl 284 IF( sn_cfctl%l_prtctl 285 IF( ln_timing 286 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 287 IF( sn_cfctl%l_prtctl ) CALL ice_prt3D ('iceupdate') ! prints 288 IF( ln_timing ) CALL timing_stop ('ice_update') ! timing 286 289 ! 287 290 END SUBROUTINE ice_update_flx -
NEMO/trunk/src/OCE/BDY/bdyice.F90
r13472 r13601 61 61 !!---------------------------------------------------------------------- 62 62 ! controls 63 IF( ln_timing ) CALL timing_start('bdy_ice_thd') ! timing 64 IF( ln_icediachk ) CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 65 IF( ln_icediachk ) CALL ice_cons2D (0,'bdy_ice_thd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 63 IF( ln_timing ) CALL timing_start('bdy_ice_thd') ! timing 66 64 ! 67 65 CALL ice_var_glo2eqv … … 110 108 ! 111 109 ! controls 112 IF( ln_icectl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints 113 IF( ln_icediachk ) CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 114 IF( ln_icediachk ) CALL ice_cons2D (1,'bdy_ice_thd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 115 IF( ln_timing ) CALL timing_stop ('bdy_ice_thd') ! timing 110 IF( ln_icectl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints 111 IF( ln_timing ) CALL timing_stop ('bdy_ice_thd') ! timing 116 112 ! 117 113 END SUBROUTINE bdy_ice
Note: See TracChangeset
for help on using the changeset viewer.
