Changeset 14717 for NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcblk.F90

Timestamp:

2021-04-16T11:42:56+02:00 (3 years ago)

Author:

clem

Message:

4.0-HEAD: correctly handle diagnostics of mass, salt and heat budgets (see ticket #2652). And fix Pierre ticket #2642

File:

• NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcblk.F90 (modified) (7 diffs)

NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcblk.F90

@@ -390 +390 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
       REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3]
-      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcptrain, zcptsnw, zcptn ! Heat content per unit mass (J/kg)
+      !!---------------------------------------------------------------------
+      !
+      ! Heat content per unit mass (J/kg)
+      zcptrain(:,:) = (      sf(jp_tair)%fnow(:,:,1)        - rt0 ) * rcp  * tmask(:,:,1)
+      zcptsnw (:,:) = ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1)
+      zcptn   (:,:) =        pst(:,:)                               * rcp  * tmask(:,:,1)
       !
       ! local scalars ( place there for vector optimisation purposes)
@@ -541 +547 @@
       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
-         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
+         &     - zevap(:,:) * zcptn(:,:)                                          &   ! remove evap heat content at SST
          &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
-         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
+         &     * zcptrain(:,:)                                                    &
          &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
-         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
+         &     * zcptsnw(:,:)
       qns(:,:) = qns(:,:) * tmask(:,:,1)
       !
@@ -819 +825 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
       REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa
-      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2
       REAL(wp), DIMENSION(jpi,jpj)     ::   ztri
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zcptrain, zcptsnw, zcptn ! Heat content per unit mass (J/kg)
       !!---------------------------------------------------------------------
       !
@@ -827 +833 @@
       !
       zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
+      !
+      ! Heat content per unit mass (J/kg)
+      zcptrain(:,:) = (      sf(jp_tair)%fnow(:,:,1)        - rt0 ) * rcp  * tmask(:,:,1)
+      zcptsnw (:,:) = ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1)
+      zcptn   (:,:) =        sst_m(:,:)                             * rcp  * tmask(:,:,1)
       !
       zztmp = 1. / ( 1. - albo )
@@ -901 +912 @@
 
       ! --- heat flux associated with emp --- !
-      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
-         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
-         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
-         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
-      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
-         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
+      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * zcptn(:,:)         & ! evap at sst
+         &          + ( tprecip(:,:) - sprecip(:,:) )   *   zcptrain(:,:)         & ! liquid precip at Tair
+         &          +   sprecip(:,:) * ( 1._wp - zsnw ) * ( zcptsnw (:,:) - rLfus ) ! solid precip at min(Tair,Tsnow)
+      qemp_ice(:,:) =   sprecip(:,:) *           zsnw   * ( zcptsnw (:,:) - rLfus ) ! solid precip (only)
 
       ! --- total solar and non solar fluxes --- !
@@ -914 +923 @@
 
       ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
-      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
+      qprec_ice(:,:) = rhos * ( zcptsnw(:,:) - rLfus )
 
       ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
@@ -947 +956 @@
 
       IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
-         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) 
-         CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average)
-         CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) )   ! heat flux from evap (cell average)
-      ENDIF
-      IF( iom_use('hflx_rain_cea') ) THEN
-         ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) )
-         CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) )   ! heat flux from rain (cell average)
-      ENDIF
-      IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN
-          WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) ;   ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
-          ELSEWHERE                             ;   ztmp(:,:) = rcp * sst_m(:,:)    
-          ENDWHERE
-          ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus ) 
-          CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average)
-          CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
-          CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice)
-      ENDIF
-      !
+         CALL iom_put( 'evap_ao_cea'  , zevap(:,:) * ( 1._wp - at_i_b(:,:) ) * tmask(:,:,1)              )   ! ice-free oce evap (cell average)
+         CALL iom_put( 'hflx_evap_cea', zevap(:,:) * ( 1._wp - at_i_b(:,:) ) * tmask(:,:,1) * zcptn(:,:) )   ! heat flux from evap (cell average)
+      ENDIF
+      IF( iom_use('rain') .OR. iom_use('rain_ao_cea') .OR. iom_use('hflx_rain_cea') ) THEN
+         CALL iom_put( 'rain'         ,   tprecip(:,:) - sprecip(:,:)                             )          ! liquid precipitation 
+         CALL iom_put( 'rain_ao_cea'  , ( tprecip(:,:) - sprecip(:,:) ) * ( 1._wp - at_i_b(:,:) ) )          ! liquid precipitation over ocean (cell average)
+         CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )                    ! heat flux from rain (cell average)
+      ENDIF
+      IF(  iom_use('snow_ao_cea')   .OR. iom_use('snow_ai_cea')      .OR. &
+         & iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN
+         CALL iom_put( 'snow_ao_cea'     , sprecip(:,:)                            * ( 1._wp - zsnw(:,:) ) ) ! Snow over ice-free ocean  (cell average)
+         CALL iom_put( 'snow_ai_cea'     , sprecip(:,:)                            *           zsnw(:,:)   ) ! Snow over sea-ice         (cell average)
+         CALL iom_put( 'hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) )                         ! heat flux from snow (cell average)
+         CALL iom_put( 'hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
+         CALL iom_put( 'hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) *           zsnw(:,:)   ) ! heat flux from snow (over ice)
+      ENDIF
+      IF( iom_use('hflx_prec_cea') ) THEN                                                                    ! heat flux from precip (cell average)
+         CALL iom_put('hflx_prec_cea' ,    sprecip(:,:)                  * ( zcptsnw (:,:) - rLfus )  &
+            &                          + ( tprecip(:,:) - sprecip(:,:) ) *   zcptrain(:,:) )
+      ENDIF
+      !
+      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea' , SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) * tmask(:,:,1) ) ! Sublimation over sea-ice (cell average)
+!!clem IF( iom_use('hflx_subl_cea') ) CALL iom_put( 'hflx_subl_cea', SUM( a_i_b(:,:,:) * qevap_ice(:,:,:), dim=3 ) * tmask(:,:,1) ) ! Heat flux from sublimation (cell average)
+
       IF(ln_ctl) THEN
          CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)