Changeset 5837 for branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90

Timestamp:

2015-10-26T15:59:39+01:00 (9 years ago)

Author:

timgraham

Message:

Upgraded to r5518 of trunk (NEMO 3.6 stable)

File:

• branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90 (modified) (7 diffs)

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90

@@ -8 +8 @@
    !!            3.0  ! 2006-08  (G. Madec)  Surface module
    !!            3.2  ! 2009-07  (C. Talandier) emp mean s spread over erp area 
+   !!            3.6  ! 2014-11  (P. Mathiot  ) add ice shelf melting
    !!----------------------------------------------------------------------
 
@@ -19 +20 @@
    USE phycst          ! physical constants
    USE sbcrnf          ! ocean runoffs
+   USE sbcisf          ! ice shelf melting contribution
    USE sbcssr          ! SS damping terms
    USE in_out_manager  ! I/O manager
@@ -57 +59 @@
       !!                =1 global mean of emp set to zero at each nn_fsbc time step
       !!                =2 annual global mean corrected from previous year
+      !!                =3 global mean of emp set to zero at each nn_fsbc time step
+      !!                   & spread out over erp area depending its sign
       !! Note: if sea ice is embedded it is taken into account when computing the budget 
       !!----------------------------------------------------------------------
@@ -81 +85 @@
             IF( kn_fwb == 1 )   WRITE(numout,*) '          instantaneously set to zero'
             IF( kn_fwb == 2 )   WRITE(numout,*) '          adjusted from previous year budget'
-         ENDIF
-         !
-         area = glob_sum( e1e2t(:,:) )           ! interior global domain surface
-         !
-#if ! defined key_lim2 &&  ! defined key_lim3 && ! defined key_cice 
+            IF( kn_fwb == 3 )   WRITE(numout,*) '          fwf set to zero and spread out over erp area'
+         ENDIF
+         !
+         IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' )
+         IF( kn_fwb == 3 .AND. ln_isfcav    )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' )
+         !
+         area = glob_sum( e1e2t(:,:) * tmask(:,:,1))           ! interior global domain surface
+         ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes
+         ! and in case of no melt, it can generate HSSW.
+         !
+#if ! defined key_lim2 &&  ! defined key_lim3 && ! defined key_cice
          snwice_mass_b(:,:) = 0.e0               ! no sea-ice model is being used : no snow+ice mass
          snwice_mass  (:,:) = 0.e0
@@ -98 +108 @@
          !
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain
+            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain
             zcoef = z_fwf * rcp
-            emp(:,:) = emp(:,:) - z_fwf 
-            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:)  ! account for change to the heat budget due to fw correction
+            emp(:,:) = emp(:,:) - z_fwf              * tmask(:,:,1)
+            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
          ENDIF
          !
@@ -132 +142 @@
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes
             zcoef = fwfold * rcp
-            emp(:,:) = emp(:,:) + fwfold
-            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:)  ! account for change to the heat budget due to fw correction
+            emp(:,:) = emp(:,:) + fwfold             * tmask(:,:,1)
+            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
          ENDIF
          !
@@ -142 +152 @@
          ENDIF
          !
+      CASE ( 3 )                             !==  global fwf set to zero and spread out over erp area  ==!
+         !
+         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
+            ztmsk_pos(:,:) = tmask_i(:,:)                      ! Select <0 and >0 area of erp
+            WHERE( erp < 0._wp )   ztmsk_pos = 0._wp
+            ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
+            !
+            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp 
+            zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
+            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges) 
+            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) - snwice_fmass(:,:) ) ) / area
+            !            
+            IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation
+                zsurf_tospread      = zsurf_pos
+                ztmsk_tospread(:,:) = ztmsk_pos(:,:)
+            ELSE                             ! spread out over <0 erp area to increase precipitation
+                zsurf_tospread      = zsurf_neg
+                ztmsk_tospread(:,:) = ztmsk_neg(:,:)
+            ENDIF
+            !
+            zsum_fwf   = glob_sum( e1e2t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
+!!gm :  zsum_fwf   = z_fwf * area   ???  it is right?  I think so....
+            z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall )
+            !                                                  ! weight to respect erp field 2D structure 
+            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
+            z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
+            !                                                  ! final correction term to apply
+            zerp_cor(:,:) = -1. * z_fwf_nsrf * zsurf_tospread * z_wgt(:,:)
+            !
+!!gm   ===>>>>  lbc_lnk should be useless as all the computation is done over the whole domain !
+            CALL lbc_lnk( zerp_cor, 'T', 1. )
+            !
+            emp(:,:) = emp(:,:) + zerp_cor(:,:)
+            qns(:,:) = qns(:,:) - zerp_cor(:,:) * rcp * sst_m(:,:)  ! account for change to the heat budget due to fw correction
+            erp(:,:) = erp(:,:) + zerp_cor(:,:)
+            !
+            IF( nprint == 1 .AND. lwp ) THEN                   ! control print
+               IF( z_fwf < 0._wp ) THEN
+                  WRITE(numout,*)'   z_fwf < 0'
+                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
+               ELSE
+                  WRITE(numout,*)'   z_fwf >= 0'
+                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
+               ENDIF
+               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv'
+               WRITE(numout,*)'   z_fwf         = ', z_fwf      ,' Kg/m2/s'
+               WRITE(numout,*)'   z_fwf_nsrf    = ', z_fwf_nsrf ,' Kg/m2/s'
+               WRITE(numout,*)'   MIN(zerp_cor) = ', MINVAL(zerp_cor) 
+               WRITE(numout,*)'   MAX(zerp_cor) = ', MAXVAL(zerp_cor) 
+            ENDIF
+         ENDIF
+         !
       CASE DEFAULT                           !==  you should never be there  ==!
-         CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1 or 2' )
+         CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1, 2 or 3' )
          !
       END SELECT