Changeset 13601


Ignore:
Timestamp:
2020-10-14T17:59:34+02:00 (3 months ago)
Author:
clem
Message:

trunk: rewrite heat budget of sea ice to make it perfectly conservative by construction. Also, activating ln_icediachk now gives an ascii file (icedrift_diagnostics.ascii) containing mass, salt and heat global conservation issues (if any). In addition, 2D drift files can be outputed (icedrift_mass…) and field_def is changed accordingly. Note that advection schemes are not yet commited since there is still a restartability issue that I do not understand.

Location:
NEMO/trunk
Files:
11 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/cfgs/SHARED/field_def_nemo-ice.xml

    r13472 r13601  
    7878          <field id="isig2"        long_name="2nd principal stress component for EVP rhg"                                                                        unit=""     /> 
    7979          <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=""     /> 
    8083          <field id="normstr"      long_name="Average normal stress in sea ice"                        standard_name="average_normal_stress"                     unit="N/m"  /> 
    8184          <field id="sheastr"      long_name="Maximum shear stress in sea ice"                         standard_name="maximum_shear_stress"                      unit="N/m"  /> 
     
    165168 
    166169          <!-- 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"    /> 
    168175 
    169176     <!-- sbcssm variables --> 
     
    459466          <field field_ref="vfxthin"          name="vfxthin" /> 
    460467    
    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"  /> 
    465468        </field_group> 
    466469 
  • NEMO/trunk/src/ICE/ice.F90

    r13472 r13601  
    392392   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw         !: snw volume variation   [m/s]  
    393393   ! 
     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   ! 
    394398   !!---------------------------------------------------------------------- 
    395399   !! * Ice conservation 
     
    495499      ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj),   &  
    496500         &      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) ) 
    498503 
    499504      ! * Ice conservation 
  • NEMO/trunk/src/ICE/icecor.F90

    r13497 r13601  
    9595               zsal = sv_i(ji,jj,jl) 
    9696               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 
    9899            END_2D 
    99100         END DO 
    100101      ENDIF 
    101       !                             !----------------------------------------------------- 
    102       CALL ice_var_zapsmall         !  Zap small values                                  ! 
    103       !                             !----------------------------------------------------- 
    104102 
     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 
    105109      !                             !----------------------------------------------------- 
    106110      IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values ! 
  • NEMO/trunk/src/ICE/icectl.F90

    r13472 r13601  
    4343   PUBLIC   ice_prt 
    4444   PUBLIC   ice_prt3D 
     45   PUBLIC   ice_drift_wri 
     46   PUBLIC   ice_drift_init 
    4547 
    4648   ! thresold rates for conservation 
     
    4951   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) 
    5052   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  
    5159    
    5260   !! * Substitutions 
     
    132140 
    133141         ! -- 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 ) 
    136144 
    137145         ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     
    156164               &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_amax 
    157165            ! 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 
    163170         ENDIF 
    164171         ! 
     
    186193      ! water flux 
    187194      ! -- 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 ) 
    189197 
    190198      ! -- 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 ) 
    192200 
    193201      ! -- 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 ) 
    197207 
    198208      ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     
    204214         IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 
    205215            &                   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 
    207218      ENDIF 
    208219      ! 
     
    725736       
    726737   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 
    727867       
    728868#else 
  • NEMO/trunk/src/ICE/icedyn_rhg_evp.F90

    r13550 r13601  
    170170      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice 
    171171      !! --- 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          
    173174      !! --- SIMIP diags 
    174175      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
     
    726727            &   ) * r1_e1e2t(ji,jj) 
    727728         zdt2 = zdt * zdt 
     729 
     730         zten_i(ji,jj) = zdt 
    728731          
    729732         ! shear**2 at T points (doc eq. A16) 
     
    741744          
    742745         ! delta at T points 
    743          zdelta(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, -zdelta(ji,jj) ) ) ! 0 if delta=0 
    745          pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch 
     746         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 
    746749 
    747750      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 ) 
    749753       
    750754      ! --- 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 ) 
    752755      pstress1_i (:,:) = zs1 (:,:) 
    753756      pstress2_i (:,:) = zs2 (:,:) 
     
    778781      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    779782 
    780       ! --- stress tensor --- ! 
    781       IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
    782          ! 
    783          ALLOCATE( zsig1(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) ) 
    784787         !          
    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          
    813846      ENDIF 
    814847 
     
    946979         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 
    947980         ! close file 
    948          IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid) 
     981         IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid) 
    949982      ENDIF 
    950983       
  • NEMO/trunk/src/ICE/icestp.F90

    r13558 r13601  
    197197         IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs  
    198198         ! 
     199         IF( ln_icediachk )             CALL ice_drift_wri( kt )      ! -- Diagnostics outputs for conservation  
     200         ! 
    199201                                        CALL ice_wri( kt )            ! -- Ice outputs  
    200202         ! 
     
    281283      ! 
    282284      CALL ice_dia_init                ! initialization for diags 
     285      ! 
     286      CALL ice_drift_init              ! initialization for diags of conservation 
    283287      ! 
    284288      fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction 
     
    341345      ENDIF 
    342346      ! 
    343       IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY') 
    344       ! 
    345347      rDt_ice   = REAL(nn_fsbc) * rn_Dt          !--- sea-ice timestep and its inverse 
    346348      r1_Dt_ice = 1._wp / rDt_ice 
     
    438440      diag_trp_ei(:,:) = 0._wp   ;   diag_trp_es(:,:) = 0._wp 
    439441      diag_trp_sv(:,:) = 0._wp 
     442      diag_adv_mass(:,:) = 0._wp 
     443      diag_adv_salt(:,:) = 0._wp 
     444      diag_adv_heat(:,:) = 0._wp 
    440445       
    441446   END SUBROUTINE diag_set0 
  • NEMO/trunk/src/ICE/icethd.F90

    r13547 r13601  
    1818   USE ice            ! sea-ice: variables 
    1919!!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_cpl 
     20   USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl 
    2121   USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 
    2222      &                 qml_ice, qcn_ice, qtr_ice_top 
     
    5252   LOGICAL ::   ln_icedO         ! activate ice growth in open-water (T) or not (F) 
    5353   LOGICAL ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F) 
    54    LOGICAL ::   ln_leadhfx       !  heat in the leads is used to melt sea-ice before warming the ocean 
     54   LOGICAL ::   ln_leadhfx       ! heat in the leads is used to melt sea-ice before warming the ocean 
    5555 
    5656   !! for convergence tests 
     
    9191      ! 
    9292      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           ! 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   ! 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) 
    9797      ! 
    9898      !!------------------------------------------------------------------- 
     
    124124               &                    (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   & 
    125125               &                     + 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) ) ) 
    126128         END_2D 
    127129      ELSE      !  if no ice dynamics => transmit directly the atmospheric stress to the ocean 
     
    130132               &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
    131133               &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 
     134            zvel(ji,jj) = 0._wp 
    132135         END_2D 
    133136      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 ) 
    135138      ! 
    136139      !--------------------------------------------------------------------! 
     
    140143         rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 
    141144         ! 
    142          !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
    143          !           !  practically no "direct lateral ablation" 
    144          !            
    145          !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean  
    146          !           !  temperature and turbulent mixing (McPhee, 1992) 
    147          ! 
    148145         ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 
    149146         zqld =  tmask(ji,jj,1) * rDt_ice *  & 
     
    151148            &      ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 
    152149 
    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) 
    154152         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) 
    155153         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) 
    158158         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          
    162161         ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach  
    163162         !                              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) ) 
    175178            qlead(ji,jj) = 0._wp 
    176179         ELSE 
    177180            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 ) 
    178188         ENDIF 
    179189         ! 
    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         ! 
    183202      END_2D 
    184203       
     
    191210      ENDIF 
    192211 
    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 - qsensible 
    197       !     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) + solar 
    199       qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean                
    200          &             - qlead(:,:) * r1_Dt_ice                                &  ! heat flux taken from the ocean where there is open water ice formation 
    201          &             - at_i (:,:) * qsb_ice_bot(:,:)                         &  ! heat flux taken by sensible flux 
    202          &             - at_i (:,:) * fhld       (:,:)                            ! heat flux taken during bottom growth/melt  
    203       !                                                                           !    (fhld should be 0 while bott growth) 
    204212      !-------------------------------------------------------------------------------------------! 
    205213      ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 
     
    418426         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res       ) 
    419427         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     ) 
    421428         ! 
    422429         ! ocean surface fields 
     
    510517         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res     ) 
    511518         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   ) 
    513519         ! 
    514520         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  
    556556         !     
    557557         ! 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_ice 
     558         !!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 
    559559 
    560560         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  
    131131 
    132132      ! 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 
    135137      END WHERE 
    136138 
     
    146148         ! 
    147149         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 landfast 
     150            IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 
    149151               ! -- Wind stress -- ! 
    150152               ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   & 
     
    198200      ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice 
    199201      !------------------------------------------------------------------------------! 
    200       ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice 
     202      ! it occurs if cooling 
    201203 
    202204      ! Identify grid points where new ice forms 
    203205      npti = 0   ;   nptidx(:) = 0 
    204206      DO_2D( 1, 1, 1, 1 ) 
    205          IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN 
     207         IF ( qlead(ji,jj)  <  0._wp ) THEN 
    206208            npti = npti + 1 
    207209            nptidx( npti ) = (jj - 1) * jpi + ji 
  • NEMO/trunk/src/ICE/iceupdate.F90

    r13497 r13601  
    2424   USE traqsr         ! add penetration of solar flux in the calculation of heat budget 
    2525   USE icectl         ! sea-ice: control prints 
    26    USE bdy_oce , ONLY : ln_bdy 
    2726   USE zdfdrg  , ONLY : ln_drgice_imp 
    2827   ! 
     
    9291      ! 
    9392      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
    94       REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
    9593      REAL(wp) ::   zqsr             ! New solar flux received by the ocean 
    9694      REAL(wp), DIMENSION(jpi,jpj) ::   z2d                  ! 2D workspace 
     
    103101         WRITE(numout,*)'~~~~~~~~~~~~~~' 
    104102      ENDIF 
     103 
     104      ! Net heat flux on top of the ice-ocean (W.m-2) 
     105      !---------------------------------------------- 
     106      qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:)  
    105107 
    106108      ! --- case we bypass ice thermodynamics --- ! 
     
    115117      DO_2D( 1, 1, 1, 1 ) 
    116118 
    117          ! Solar heat flux reaching the ocean = zqsr (W.m-2)  
     119         ! Solar heat flux reaching the ocean (max) = zqsr (W.m-2)  
    118120         !--------------------------------------------------- 
    119121         zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) ) 
     
    121123         ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2)  
    122124         !--------------------------------------------------- 
    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          
    131130         ! New qsr and qns used to compute the oceanic heat flux at the next time step 
    132131         !---------------------------------------------------------------------------- 
    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 
    134146         qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr               
    135  
     147          
    136148         ! Mass flux at the atm. surface        
    137149         !----------------------------------- 
     
    140152         ! Mass flux at the ocean surface       
    141153         !------------------------------------ 
    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 
    151155         wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   & 
    152156            &           + 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 ocean 
     157          
     158         ! snw-ocean mass flux 
    155159         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 
    161164 
    162165         ! Salt flux at the ocean surface       
     
    262265      CALL iom_put ('hfxdif'     , hfx_dif     )   ! heat flux used for ice temperature change 
    263266      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 
    265268 
    266269      ! heat fluxes associated with mass exchange (freeze/melt/precip...) 
     
    279282      !--------- 
    280283#if ! defined key_agrif 
    281       IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation 
     284      IF( ln_icediachk      )   CALL ice_cons_final('iceupdate')                                       ! conservation 
    282285#endif 
    283       IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 
    284       IF( sn_cfctl%l_prtctl              )   CALL ice_prt3D     ('iceupdate')                                       ! prints 
    285       IF( ln_timing                      )   CALL timing_stop   ('ice_update')                                      ! 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 
    286289      ! 
    287290   END SUBROUTINE ice_update_flx 
  • NEMO/trunk/src/OCE/BDY/bdyice.F90

    r13472 r13601  
    6161      !!---------------------------------------------------------------------- 
    6262      ! 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 
    6664      ! 
    6765      CALL ice_var_glo2eqv 
     
    110108      ! 
    111109      ! 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 
    116112      ! 
    117113   END SUBROUTINE bdy_ice 
Note: See TracChangeset for help on using the changeset viewer.