Changeset 13589


Ignore:
Timestamp:
2020-10-14T15:35:49+02:00 (3 months ago)
Author:
clem
Message:

4.0-HEAD: 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…) but since I did not want to change the xml files, one has to uncomment some lines in icectl.F90 and add the corresponding fields in field_def_oce.xml

Location:
NEMO/releases/r4.0/r4.0-HEAD/src
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/ice.F90

    r13444 r13589  
    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/releases/r4.0/r4.0-HEAD/src/ICE/icecor.F90

    r13284 r13589  
    9696                  zsal = sv_i(ji,jj,jl) 
    9797                  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)  ) 
    98                   sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux 
     98                  IF( kn /= 0 ) & ! no ice-ocean exchanges if kn=0 (for bdy for instance) otherwise conservation diags will fail 
     99                     &   sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux 
    99100               END DO 
    100101            END DO 
    101102         END DO 
    102103      ENDIF 
    103       !                             !----------------------------------------------------- 
    104       CALL ice_var_zapsmall         !  Zap small values                                  ! 
    105       !                             !----------------------------------------------------- 
    106104 
     105      IF( kn /= 0 ) THEN   ! no zapsmall if kn=0 (for bdy for instance) because we do not want ice-ocean exchanges (wfx,sfx,hfx) 
     106         !                                                              otherwise conservation diags will fail 
     107         !                          !----------------------------------------------------- 
     108         CALL ice_var_zapsmall      !  Zap small values                                  ! 
     109         !                          !----------------------------------------------------- 
     110      ENDIF 
    107111      !                             !----------------------------------------------------- 
    108112      IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values ! 
     
    117121            END DO 
    118122         END DO 
    119          CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1., v_ice, 'V', -1. ) 
     123         CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1._wp, v_ice, 'V', -1._wp ) 
    120124      ENDIF 
    121125 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icectl.F90

    r13284 r13589  
    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      ! 
     
    747758       
    748759   END SUBROUTINE ice_prt3D 
    749        
     760 
     761 
     762   SUBROUTINE ice_drift_wri( kt ) 
     763      !!------------------------------------------------------------------- 
     764      !!                     ***  ROUTINE ice_drift_wri *** 
     765      !! 
     766      !! ** Purpose : conservation of mass, salt and heat 
     767      !!              write the drift in a ascii file at each time step 
     768      !!              and the total run drifts 
     769      !!------------------------------------------------------------------- 
     770      INTEGER, INTENT(in) ::   kt   ! ice time-step index 
     771      ! 
     772      INTEGER  ::   ji, jj 
     773      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 
     774      !!REAL(wp), DIMENSION(jpi,jpj) ::   zdiag_mass2D, zdiag_salt2D, zdiag_heat2D 
     775      !!------------------------------------------------------------------- 
     776      ! 
     777      IF( kt == nit000 .AND. lwp ) THEN 
     778         WRITE(numout,*) 
     779         WRITE(numout,*) 'ice_drift_wri: sea-ice drifts' 
     780         WRITE(numout,*) '~~~~~~~~~~~~~' 
     781      ENDIF 
     782      ! 
     783      !clem: the following lines check the ice drift in 2D. 
     784      !      to use this check, uncomment those lines and add the 3 fields in field_def_ice.xml 
     785      !!! 2D budgets (must be close to 0) 
     786      !!IF( iom_use('icedrift_mass') .OR. iom_use('icedrift_salt') .OR. iom_use('icedrift_heat') ) THEN 
     787      !!   DO jj = 1, jpj 
     788      !!      DO ji = 1, jpi 
     789      !!         zdiag_mass2D(ji,jj) =   wfx_ice(ji,jj)   + wfx_snw(ji,jj)   + wfx_spr(ji,jj) + wfx_sub(ji,jj) & 
     790      !!            &                  + diag_vice(ji,jj) + diag_vsnw(ji,jj) - diag_adv_mass(ji,jj) 
     791      !!         zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) 
     792      !!         zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj) 
     793      !!      END DO 
     794      !!   END DO 
     795      !!   ! 
     796      !!   ! write outputs 
     797      !!   CALL iom_put( 'icedrift_mass', zdiag_mass2D ) 
     798      !!   CALL iom_put( 'icedrift_salt', zdiag_salt2D ) 
     799      !!   CALL iom_put( 'icedrift_heat', zdiag_heat2D ) 
     800      !!ENDIF 
     801 
     802      ! -- mass diag -- ! 
     803      zdiag_mass     = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub & 
     804         &                                  + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rdt_ice 
     805      zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rdt_ice 
     806 
     807      ! -- salt diag -- ! 
     808      zdiag_salt     = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rdt_ice * 1.e-3 
     809      zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rdt_ice * 1.e-3 
     810 
     811      ! -- heat diag -- ! 
     812      zdiag_heat     = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 
     813      zdiag_adv_heat = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 
     814 
     815      !                    ! write out to file 
     816      IF( lwp ) THEN 
     817         ! check global drift (must be close to 0) 
     818         WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift     [kg]', zdiag_mass 
     819         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift     [kg]', zdiag_salt 
     820         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift     [W] ', zdiag_heat 
     821         ! check drift from advection scheme (can be /=0 with bdy but not sure why) 
     822         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'mass drift adv [kg]', zdiag_adv_mass 
     823         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift adv [kg]', zdiag_adv_salt 
     824         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift adv [W] ', zdiag_adv_heat 
     825      ENDIF 
     826      !                    ! drifts 
     827      rdiag_icemass = rdiag_icemass + zdiag_mass 
     828      rdiag_icesalt = rdiag_icesalt + zdiag_salt 
     829      rdiag_iceheat = rdiag_iceheat + zdiag_heat 
     830      rdiag_adv_icemass = rdiag_adv_icemass + zdiag_adv_mass 
     831      rdiag_adv_icesalt = rdiag_adv_icesalt + zdiag_adv_salt 
     832      rdiag_adv_iceheat = rdiag_adv_iceheat + zdiag_adv_heat 
     833      ! 
     834      !                    ! output drifts and close ascii file 
     835      IF( kt == nitend - nn_fsbc + 1 .AND. lwp ) THEN 
     836         ! to ascii file 
     837         WRITE(numicedrift,*) '******************************************' 
     838         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift     [kg]', rdiag_icemass 
     839         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift adv [kg]', rdiag_adv_icemass 
     840         WRITE(numicedrift,*) '******************************************' 
     841         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift     [kg]', rdiag_icesalt 
     842         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift adv [kg]', rdiag_adv_icesalt 
     843         WRITE(numicedrift,*) '******************************************' 
     844         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift     [W] ', rdiag_iceheat 
     845         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift adv [W] ', rdiag_adv_iceheat 
     846         CLOSE( numicedrift ) 
     847         ! 
     848         ! to ocean output 
     849         WRITE(numout,*) 
     850         WRITE(numout,*) 'ice_drift_wri: ice drifts information for the run ' 
     851         WRITE(numout,*) '~~~~~~~~~~~~~' 
     852         ! check global drift (must be close to 0) 
     853         WRITE(numout,*) '   sea-ice mass drift     [kg] = ', rdiag_icemass 
     854         WRITE(numout,*) '   sea-ice salt drift     [kg] = ', rdiag_icesalt 
     855         WRITE(numout,*) '   sea-ice heat drift     [W]  = ', rdiag_iceheat 
     856         ! check drift from advection scheme (can be /=0 with bdy but not sure why) 
     857         WRITE(numout,*) '   sea-ice mass drift adv [kg] = ', rdiag_adv_icemass 
     858         WRITE(numout,*) '   sea-ice salt drift adv [kg] = ', rdiag_adv_icesalt 
     859         WRITE(numout,*) '   sea-ice heat drift adv [W]  = ', rdiag_adv_iceheat 
     860      ENDIF 
     861      ! 
     862   END SUBROUTINE ice_drift_wri 
     863 
     864   SUBROUTINE ice_drift_init 
     865      !!---------------------------------------------------------------------- 
     866      !!                  ***  ROUTINE ice_drift_init  *** 
     867      !!                    
     868      !! ** Purpose :   create output file, initialise arrays 
     869      !!---------------------------------------------------------------------- 
     870      ! 
     871      IF( .NOT.ln_icediachk ) RETURN ! exit 
     872      ! 
     873      IF(lwp) THEN 
     874         WRITE(numout,*) 
     875         WRITE(numout,*) 'ice_drift_init: Output ice drifts to ',TRIM(clname), ' file' 
     876         WRITE(numout,*) '~~~~~~~~~~~~~' 
     877         WRITE(numout,*) 
     878         ! 
     879         ! create output ascii file 
     880         CALL ctl_opn( numicedrift, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, narea ) 
     881         WRITE(numicedrift,*) 'Timestep  Drifts' 
     882         WRITE(numicedrift,*) '******************************************' 
     883      ENDIF 
     884      ! 
     885      rdiag_icemass = 0._wp 
     886      rdiag_icesalt = 0._wp 
     887      rdiag_iceheat = 0._wp 
     888      rdiag_adv_icemass = 0._wp 
     889      rdiag_adv_icesalt = 0._wp 
     890      rdiag_adv_iceheat = 0._wp 
     891      ! 
     892   END SUBROUTINE ice_drift_init 
     893 
    750894#else 
    751895   !!---------------------------------------------------------------------- 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_pra.F90

    r13564 r13589  
    8888      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    8989      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    90       REAL(wp) ::   zdt                     !   -      - 
     90      REAL(wp) ::   zdt, z1_dt              !   -      - 
    9191      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
    9292      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
     
    100100      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es 
    101101      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei 
     102      !! diagnostics 
     103      REAL(wp), DIMENSION(jpi,jpj)            ::   zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat       
    102104      !!---------------------------------------------------------------------- 
    103105      ! 
     
    185187      ENDIF 
    186188      zdt = rdt_ice / REAL(icycle) 
     189      z1_dt = 1._wp / zdt 
    187190       
    188191      ! --- transport --- ! 
     
    191194 
    192195      DO jt = 1, icycle 
     196 
     197         ! diagnostics 
     198         zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     199         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
     200         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     201            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 
    193202 
    194203         ! record at_i before advection (for open water) 
     
    351360         END DO 
    352361         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1. ) 
     362         ! 
     363         ! --- diagnostics --- ! 
     364         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     365            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
     366         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
     367            &                                        - zdiag_adv_salt(:,:) ) * z1_dt 
     368         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     369            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 
     370            &                                        - zdiag_adv_heat(:,:) ) * z1_dt 
    353371         ! 
    354372         ! --- Ensure non-negative fields --- ! 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_umx.F90

    r13566 r13589  
    9292      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    9393      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    94       REAL(wp) ::   zdt, zvi_cen 
     94      REAL(wp) ::   zdt, z1_dt, zvi_cen 
    9595      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
    9696      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx, zcu_box, zcv_box 
     
    104104      ! 
    105105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs  
     106      !! diagnostics 
     107      REAL(wp), DIMENSION(jpi,jpj)            ::   zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat       
    106108      !!---------------------------------------------------------------------- 
    107109      ! 
     
    189191      ENDIF 
    190192      zdt = rdt_ice / REAL(icycle) 
    191  
     193      z1_dt = 1._wp / zdt 
     194       
    192195      ! --- transport --- ! 
    193196      zudy(:,:) = pu_ice(:,:) * e2u(:,:) 
     
    219222      !---------------! 
    220223      DO jt = 1, icycle 
     224 
     225         ! diagnostics 
     226         zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     227         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
     228         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     229            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 
    221230 
    222231         ! record at_i before advection (for open water) 
     
    414423         END DO 
    415424         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp ) 
     425         ! 
     426         ! --- diagnostics --- ! 
     427         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     428            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
     429         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
     430            &                                        - zdiag_adv_salt(:,:) ) * z1_dt 
     431         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     432            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 
     433            &                                        - zdiag_adv_heat(:,:) ) * z1_dt 
    416434         ! 
    417435         ! --- Ensure non-negative fields and in-bound thicknesses --- ! 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_rhg_evp.F90

    r13549 r13589  
    990990         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 
    991991         ! close file 
    992          IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid) 
     992         IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid) 
    993993      ENDIF 
    994994       
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icestp.F90

    r13445 r13589  
    198198         IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs  
    199199         ! 
     200         IF( ln_icediachk )             CALL ice_drift_wri( kt )      ! -- Diagnostics outputs for conservation  
     201         ! 
    200202                                        CALL ice_wri( kt )            ! -- Ice outputs  
    201203         ! 
     
    252254      END WHERE 
    253255      ! 
     256      CALL diag_set0                   ! set diag of mass, heat and salt fluxes to 0: needed for Agrif child grids 
     257      ! 
    254258      CALL ice_itd_init                ! ice thickness distribution initialization 
    255259      ! 
    256260      CALL ice_thd_init                ! set ice thermodynics parameters (clem: important to call it first for melt ponds) 
    257261      ! 
    258       !                                ! Initial sea-ice state 
    259       CALL ice_istate_init 
     262      CALL ice_sbc_init                ! set ice-ocean and ice-atm. coupling parameters 
     263      ! 
     264      CALL ice_istate_init             ! Initial sea-ice state 
    260265      IF ( ln_rstart .OR. nn_iceini_file == 2 ) THEN 
    261266         CALL ice_rst_read                      ! start from a restart file 
     
    266271      CALL ice_var_agg(1) 
    267272      ! 
    268       CALL ice_sbc_init                ! set ice-ocean and ice-atm. coupling parameters 
    269       ! 
    270273      CALL ice_dyn_init                ! set ice dynamics parameters 
    271274      ! 
     
    275278      ! 
    276279      CALL ice_dia_init                ! initialization for diags 
     280      ! 
     281      CALL ice_drift_init              ! initialization for diags of conservation 
    277282      ! 
    278283      fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction 
     
    337342      ENDIF 
    338343      ! 
    339       IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY') 
    340       ! 
    341344      rdt_ice   = REAL(nn_fsbc) * rdt          !--- sea-ice timestep and its inverse 
    342345      r1_rdtice = 1._wp / rdt_ice 
     
    434437      diag_trp_ei(:,:) = 0._wp   ;   diag_trp_es(:,:) = 0._wp 
    435438      diag_trp_sv(:,:) = 0._wp 
     439      diag_adv_mass(:,:) = 0._wp 
     440      diag_adv_salt(:,:) = 0._wp 
     441      diag_adv_heat(:,:) = 0._wp 
    436442       
    437443   END SUBROUTINE diag_set0 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd.F90

    r13284 r13589  
    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      !!------------------------------------------------------------------- 
     
    125125                  &                    (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   & 
    126126                  &                     + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1) 
     127               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) ) + & 
     128                  &                         ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 
    127129            END DO 
    128130         END DO 
    129       ELSE      !  if no ice dynamics => transmit directly the atmospheric stress to the ocean 
     131      ELSE      !  if no ice dynamics => transfer directly the atmospheric stress to the ocean 
    130132         DO jj = 2, jpjm1 
    131133            DO ji = fs_2, fs_jpim1 
     
    133135                  &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
    134136                  &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 
     137               zvel(ji,jj) = 0._wp 
    135138            END DO 
    136139         END DO 
    137140      ENDIF 
    138       CALL lbc_lnk( 'icethd', zfric, 'T',  1. ) 
     141      CALL lbc_lnk_multi( 'icethd', zfric, 'T',  1._wp, zvel, 'T', 1._wp ) 
    139142      ! 
    140143      !--------------------------------------------------------------------! 
     
    145148            rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 
    146149            ! 
    147             !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
    148             !           !  practically no "direct lateral ablation" 
    149             !            
    150             !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean  
    151             !           !  temperature and turbulent mixing (McPhee, 1992) 
    152             ! 
    153150            ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 
    154151            zqld =  tmask(ji,jj,1) * rdt_ice *  & 
     
    156153               &      ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 
    157154 
    158             ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- ! 
     155            ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 
     156            !     (mostly<0 but >0 if supercooling) 
    159157            zqfr     = rau0 * 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) 
    160158            zqfr_neg = MIN( zqfr , 0._wp )                                                                    ! only < 0 
    161  
    162             ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2) 
     159            zqfr_pos = MAX( zqfr , 0._wp )                                                                    ! only > 0 
     160 
     161            ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 
     162            !     (mostly>0 but <0 if supercooling) 
    163163            zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin )  
    164             qsb_ice_bot(ji,jj) = rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 
    165  
    166             qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 
     164            qsb_ice_bot(ji,jj) = rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 
     165 
    167166            ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach  
    168167            !                              the freezing point, so that we do not have SST < T_freeze 
    169             !                              This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0 
    170  
    171             !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice 
    172             qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr ) 
    173  
    174             ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting  
    175             ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 
    176             IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 
    177                IF( ln_leadhfx ) THEN   ;   fhld(ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
    178                ELSE                    ;   fhld(ji,jj) = 0._wp 
    179                ENDIF 
     168            !                              This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 
     169            !                              The following formulation is ok for both normal conditions and supercooling 
     170            qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 
     171 
     172            ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 
     173            !     qlead is the energy received from the atm. in the leads. 
     174            !     If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld  (W/m2) 
     175            !     If cooling (zqld <  0), then the energy in the leads is used to grow ice in open water    => qlead (J.m-2) 
     176            IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
     177               ! upper bound for fhld: fhld should be equal to zqld 
     178               !                        but we have to make sure that this heat will not make the sst drop below the freezing point 
     179               !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 
     180               !                        The following formulation is ok for both normal conditions and supercooling 
     181               fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) &  ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
     182                  &                                 - qsb_ice_bot(ji,jj) ) 
    180183               qlead(ji,jj) = 0._wp 
    181184            ELSE 
    182185               fhld (ji,jj) = 0._wp 
     186               ! upper bound for qlead: qlead should be equal to zqld 
     187               !                        but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 
     188               !                        The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 
     189               !                        and freezing point is reached if zqfr = zqld - qsb*a/dt 
     190               !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 
     191               !                        The following formulation is ok for both normal conditions and supercooling 
     192               qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr ) 
    183193            ENDIF 
    184194            ! 
    185             ! Net heat flux on top of the ice-ocean [W.m-2] 
    186             ! --------------------------------------------- 
    187             qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)  
     195            ! If ice is landfast and ice concentration reaches its max 
     196            ! => stop ice formation in open water 
     197            IF(  zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 )   qlead(ji,jj) = 0._wp 
     198            ! 
     199            ! If the grid cell is almost fully covered by ice (no leads) 
     200            ! => stop ice formation in open water 
     201            IF( at_i(ji,jj) >= (1._wp - epsi10) )   qlead(ji,jj) = 0._wp 
     202            ! 
     203            ! If ln_leadhfx is false 
     204            ! => do not use energy of the leads to melt sea-ice 
     205            IF( .NOT.ln_leadhfx )   fhld(ji,jj) = 0._wp 
     206            ! 
    188207         END DO 
    189208      END DO 
     
    197216      ENDIF 
    198217 
    199       ! --------------------------------------------------------------------- 
    200       ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2] 
    201       ! --------------------------------------------------------------------- 
    202       !     First  step here              :  non solar + precip - qlead - qsensible 
    203       !     Second step in icethd_dh      :  heat remaining if total melt (zq_rema)  
    204       !     Third  step in iceupdate.F90  :  heat from ice-ocean mass exchange (zf_mass) + solar 
    205       qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean                
    206          &             - qlead(:,:) * r1_rdtice                                &  ! heat flux taken from the ocean where there is open water ice formation 
    207          &             - at_i (:,:) * qsb_ice_bot(:,:)                         &  ! heat flux taken by sensible flux 
    208          &             - at_i (:,:) * fhld       (:,:)                            ! heat flux taken during bottom growth/melt  
    209       !                                                                           !    (fhld should be 0 while bott growth) 
    210218      !-------------------------------------------------------------------------------------------! 
    211219      ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 
     
    425433         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res       ) 
    426434         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif   ) 
    427          CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai     ) 
    428435         ! 
    429436         ! ocean surface fields 
     
    517524         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res     ) 
    518525         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 
    519          CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai   ) 
    520526         ! 
    521527         CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d    (1:npti), qns_ice    (:,:,kl) ) 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd_dh.F90

    r13284 r13589  
    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_rdtice 
     558         !!!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 
    559559 
    560560         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd_do.F90

    r11536 r13589  
    129129 
    130130      ! Default new ice thickness 
    131       WHERE( qlead(:,:) < 0._wp  .AND. tau_icebfr(:,:) == 0._wp )   ;   ht_i_new(:,:) = rn_hinew ! if cooling and no landfast 
    132       ELSEWHERE                                                     ;   ht_i_new(:,:) = 0._wp 
     131      WHERE( qlead(:,:) < 0._wp ) ! cooling 
     132         ht_i_new(:,:) = rn_hinew 
     133      ELSEWHERE 
     134         ht_i_new(:,:) = 0._wp 
    133135      END WHERE 
    134136 
     
    145147         DO jj = 2, jpjm1 
    146148            DO ji = 2, jpim1 
    147                IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast 
     149               IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 
    148150                  ! -- Wind stress -- ! 
    149151                  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 
     
    204206      DO jj = 1, jpj 
    205207         DO ji = 1, jpi 
    206             IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN 
     208            IF ( qlead(ji,jj) < 0._wp ) THEN 
    207209               npti = npti + 1 
    208210               nptidx( npti ) = (jj - 1) * jpi + ji 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceupdate.F90

    r13284 r13589  
    2525   USE traqsr         ! add penetration of solar flux in the calculation of heat budget 
    2626   USE icectl         ! sea-ice: control prints 
    27    USE bdy_oce , ONLY : ln_bdy 
    2827   USE zdfdrg  , ONLY : ln_drgice_imp 
    2928   ! 
     
    9392      ! 
    9493      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
    95       REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
    9694      REAL(wp) ::   zqsr             ! New solar flux received by the ocean 
    9795      REAL(wp), DIMENSION(jpi,jpj) ::   z2d                  ! 2D workspace 
     
    104102         WRITE(numout,*)'~~~~~~~~~~~~~~' 
    105103      ENDIF 
     104 
     105      ! Net heat flux on top of the ice-ocean (W.m-2) 
     106      !---------------------------------------------- 
     107      qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:)  
    106108 
    107109      ! --- case we bypass ice thermodynamics --- ! 
     
    117119         DO ji = 1, jpi 
    118120 
    119             ! Solar heat flux reaching the ocean = zqsr (W.m-2)  
     121            ! Solar heat flux reaching the ocean (max) = zqsr (W.m-2)  
    120122            !--------------------------------------------------- 
    121123            zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) ) 
     
    123125            ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2)  
    124126            !--------------------------------------------------- 
    125             zqmass           = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 
    126             qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + zqmass + zqsr 
    127  
    128             ! Add the residual from heat diffusion equation and sublimation (W.m-2) 
    129             !---------------------------------------------------------------------- 
    130             qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + hfx_err_dif(ji,jj) +   & 
    131                &             ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) ) 
    132  
     127            qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) & 
     128               &                                - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) & 
     129               &                                + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) & 
     130               &                                + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj)                  
     131             
    133132            ! New qsr and qns used to compute the oceanic heat flux at the next time step 
    134133            !---------------------------------------------------------------------------- 
    135             qsr(ji,jj) = zqsr                                       
     134            ! if warming and some ice remains, then we suppose that the whole solar flux has been consumed to melt the ice 
     135            ! else ( cooling or no ice left ), then we suppose that     no    solar flux has been consumed 
     136            ! 
     137            IF( fhld(ji,jj) > 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN   !-- warming and some ice remains 
     138               !                                        solar flux transmitted thru the 1st level of the ocean (i.e. not used by sea-ice) 
     139               qsr(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * ( 1._wp - frq_m(ji,jj) ) & 
     140                  !                                   + solar flux transmitted thru ice (also not used by sea-ice) 
     141                  &             + SUM( a_i_b(ji,jj,:) * qtr_ice_bot(ji,jj,:) ) 
     142               ! 
     143            ELSE                                                       !-- cooling or no ice left 
     144               qsr(ji,jj) = zqsr 
     145            ENDIF 
     146            ! 
     147            ! the non-solar is simply derived from the solar flux 
    136148            qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr               
    137149 
     
    142154            ! Mass flux at the ocean surface       
    143155            !------------------------------------ 
    144             !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED) 
    145             !  -------------------------------------------------------------------------------------  
    146             !  The idea of this approach is that the system that we consider is the ICE-OCEAN system 
    147             !  Thus  FW  flux  =  External ( E-P+snow melt) 
    148             !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW 
    149             !                     Associated to Ice formation AND Ice melting 
    150             !                     Even if i see Ice melting as a FW and SALT flux 
    151             !         
    152             ! mass flux from ice/ocean 
     156            ! ice-ocean  mass flux 
    153157            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   & 
    154158               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 
    155159 
    156             ! add the snow melt water to snow mass flux to the ocean 
     160            ! snw-ocean mass flux 
    157161            wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj) 
    158162 
    159             ! mass flux at the ocean/ice interface 
    160             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 
    161             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) 
    162  
     163            ! total mass flux at the ocean/ice interface 
     164            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 
     165            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 
    163166 
    164167            ! Salt flux at the ocean surface       
     
    265268      CALL iom_put ('hfxdif'     , hfx_dif     )   ! heat flux used for ice temperature change 
    266269      CALL iom_put ('hfxsnw'     , hfx_snw     )   ! heat flux used for snow melt  
    267       CALL iom_put ('hfxerr'     , hfx_err_dif )   ! heat flux error after heat diffusion (included in qt_oce_ai) 
     270      CALL iom_put ('hfxerr'     , hfx_err_dif )   ! heat flux error after heat diffusion 
    268271 
    269272      ! heat fluxes associated with mass exchange (freeze/melt/precip...) 
     
    282285      !--------- 
    283286#if ! defined key_agrif 
    284       IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation 
     287      IF( ln_icediachk )   CALL ice_cons_final('iceupdate')                                       ! conservation 
    285288#endif 
    286       IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 
    287       IF( ln_ctl                         )   CALL ice_prt3D     ('iceupdate')                                       ! prints 
    288       IF( ln_timing                      )   CALL timing_stop   ('ice_update')                                      ! timing 
     289      IF( ln_icectl    )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 
     290      IF( ln_ctl       )   CALL ice_prt3D     ('iceupdate')                                       ! prints 
     291      IF( ln_timing    )   CALL timing_stop   ('ice_update')                                      ! timing 
    289292      ! 
    290293   END SUBROUTINE ice_update_flx 
  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/BDY/bdyice.F90

    r13284 r13589  
    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.