source: NEMO/trunk/src/OCE/SBC/sbcfwb.F90 @ 10570

Last change on this file since 10570 was 10570, checked in by acc, 20 months ago

Trunk update to implement finer control over the choice of text report files generated. See ticket: #2167

  • Property svn:keywords set to Id
File size: 12.1 KB
Line 
1MODULE sbcfwb
2   !!======================================================================
3   !!                       ***  MODULE  sbcfwb  ***
4   !! Ocean fluxes   : domain averaged freshwater budget
5   !!======================================================================
6   !! History :  OPA  ! 2001-02  (E. Durand)  Original code
7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
8   !!            3.0  ! 2006-08  (G. Madec)  Surface module
9   !!            3.2  ! 2009-07  (C. Talandier) emp mean s spread over erp area
10   !!            3.6  ! 2014-11  (P. Mathiot  ) add ice shelf melting
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   sbc_fwb       : freshwater budget for global ocean configurations (free surface & forced mode)
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers
17   USE dom_oce        ! ocean space and time domain
18   USE sbc_oce        ! surface ocean boundary condition
19   USE sbc_ice , ONLY : snwice_mass, snwice_mass_b, snwice_fmass
20   USE phycst         ! physical constants
21   USE sbcrnf         ! ocean runoffs
22   USE sbcisf         ! ice shelf melting contribution
23   USE sbcssr         ! Sea-Surface damping terms
24   !
25   USE in_out_manager ! I/O manager
26   USE lib_mpp        ! distribued memory computing library
27   USE timing         ! Timing
28   USE lbclnk         ! ocean lateral boundary conditions
29   USE lib_fortran    !
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   sbc_fwb    ! routine called by step
35
36   REAL(wp) ::   a_fwb_b   ! annual domain averaged freshwater budget
37   REAL(wp) ::   a_fwb     ! for 2 year before (_b) and before year.
38   REAL(wp) ::   fwfold    ! fwfold to be suppressed
39   REAL(wp) ::   area      ! global mean ocean surface (interior domain)
40
41   !! * Substitutions
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
45   !! $Id$
46   !! Software governed by the CeCILL license (see ./LICENSE)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE sbc_fwb( kt, kn_fwb, kn_fsbc )
51      !!---------------------------------------------------------------------
52      !!                  ***  ROUTINE sbc_fwb  ***
53      !!
54      !! ** Purpose :   Control the mean sea surface drift
55      !!
56      !! ** Method  :   several ways  depending on kn_fwb
57      !!                =0 no control
58      !!                =1 global mean of emp set to zero at each nn_fsbc time step
59      !!                =2 annual global mean corrected from previous year
60      !!                =3 global mean of emp set to zero at each nn_fsbc time step
61      !!                   & spread out over erp area depending its sign
62      !! Note: if sea ice is embedded it is taken into account when computing the budget
63      !!----------------------------------------------------------------------
64      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
65      INTEGER, INTENT( in ) ::   kn_fsbc  !
66      INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index
67      !
68      INTEGER  ::   inum, ikty, iyear     ! local integers
69      REAL(wp) ::   z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp                ! local scalars
70      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread, zcoef          !   -      -
71      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztmsk_neg, ztmsk_pos, z_wgt ! 2D workspaces
72      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztmsk_tospread, zerp_cor    !   -      -
73      REAL(wp)   ,DIMENSION(1) ::   z_fwfprv 
74      COMPLEX(wp),DIMENSION(1) ::   y_fwfnow 
75      !!----------------------------------------------------------------------
76      !
77      IF( kt == nit000 ) THEN
78         IF(lwp) THEN
79            WRITE(numout,*)
80            WRITE(numout,*) 'sbc_fwb : FreshWater Budget correction'
81            WRITE(numout,*) '~~~~~~~'
82            IF( kn_fwb == 1 )   WRITE(numout,*) '          instantaneously set to zero'
83            IF( kn_fwb == 2 )   WRITE(numout,*) '          adjusted from previous year budget'
84            IF( kn_fwb == 3 )   WRITE(numout,*) '          fwf set to zero and spread out over erp area'
85         ENDIF
86         !
87         IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' )
88         IF( kn_fwb == 3 .AND. ln_isfcav    )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' )
89         !
90         area = glob_sum( 'sbcfwb', e1e2t(:,:) * tmask(:,:,1))           ! interior global domain surface
91         ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes
92         ! and in case of no melt, it can generate HSSW.
93         !
94#if ! defined key_si3 && ! defined key_cice
95         snwice_mass_b(:,:) = 0.e0               ! no sea-ice model is being used : no snow+ice mass
96         snwice_mass  (:,:) = 0.e0
97#endif
98         !
99      ENDIF
100
101      SELECT CASE ( kn_fwb )
102      !
103      CASE ( 1 )                             !==  global mean fwf set to zero  ==!
104         !
105         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
106            y_fwfnow(1) = local_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) )
107            CALL mpp_delay_sum( 'sbcfwb', 'fwb', y_fwfnow(:), z_fwfprv(:), kt == nitend - nn_fsbc + 1 )
108            z_fwfprv(1) = z_fwfprv(1) / area
109            zcoef = z_fwfprv(1) * rcp
110            emp(:,:) = emp(:,:) - z_fwfprv(1)        * tmask(:,:,1)
111            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
112         ENDIF
113         !
114      CASE ( 2 )                             !==  fwf budget adjusted from the previous year  ==!
115         !
116         IF( kt == nit000 ) THEN                      ! initialisation
117            !                                         ! Read the corrective factor on precipitations (fwfold)
118            CALL ctl_opn( inum, 'EMPave_old.dat', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
119            READ ( inum, "(24X,I8,2ES24.16)" ) iyear, a_fwb_b, a_fwb
120            CLOSE( inum )
121            fwfold = a_fwb                            ! current year freshwater budget correction
122            !                                         ! estimate from the previous year budget
123            IF(lwp)WRITE(numout,*)
124            IF(lwp)WRITE(numout,*)'sbc_fwb : year = ',iyear  , ' freshwater budget correction = ', fwfold
125            IF(lwp)WRITE(numout,*)'          year = ',iyear-1, ' freshwater budget read       = ', a_fwb
126            IF(lwp)WRITE(numout,*)'          year = ',iyear-2, ' freshwater budget read       = ', a_fwb_b
127         ENDIF   
128         !                                         ! Update fwfold if new year start
129         ikty = 365 * 86400 / rdt                  !!bug  use of 365 days leap year or 360d year !!!!!!!
130         IF( MOD( kt, ikty ) == 0 ) THEN
131            a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
132                                                      ! sum over the global domain
133            a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) )
134            a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s
135!!gm        !                                                      !!bug 365d year
136            fwfold =  a_fwb                           ! current year freshwater budget correction
137            !                                         ! estimate from the previous year budget
138         ENDIF
139         !
140         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes
141            zcoef = fwfold * rcp
142            emp(:,:) = emp(:,:) + fwfold             * tmask(:,:,1)
143            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
144         ENDIF
145         !
146         IF( kt == nitend .AND. lwm ) THEN            ! save fwfold value in a file (only one required)
147            CALL ctl_opn( inum, 'EMPave.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
148            WRITE( inum, "(24X,I8,2ES24.16)" ) nyear, a_fwb_b, a_fwb
149            CLOSE( inum )
150         ENDIF
151         !
152      CASE ( 3 )                             !==  global fwf set to zero and spread out over erp area  ==!
153         !
154         ALLOCATE( ztmsk_neg(jpi,jpj) , ztmsk_pos(jpi,jpj) , ztmsk_tospread(jpi,jpj) , z_wgt(jpi,jpj) , zerp_cor(jpi,jpj) )
155         !
156         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
157            ztmsk_pos(:,:) = tmask_i(:,:)                      ! Select <0 and >0 area of erp
158            WHERE( erp < 0._wp )   ztmsk_pos = 0._wp
159            ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
160            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges)
161            z_fwf     = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area
162            !           
163            IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation
164               zsurf_pos = glob_sum( 'sbcfwb', e1e2t(:,:)*ztmsk_pos(:,:) )
165               zsurf_tospread      = zsurf_pos
166               ztmsk_tospread(:,:) = ztmsk_pos(:,:)
167            ELSE                             ! spread out over <0 erp area to increase precipitation
168               zsurf_neg = glob_sum( 'sbcfwb', e1e2t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp
169               zsurf_tospread      = zsurf_neg
170               ztmsk_tospread(:,:) = ztmsk_neg(:,:)
171            ENDIF
172            !
173            zsum_fwf   = glob_sum( 'sbcfwb', e1e2t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
174!!gm :  zsum_fwf   = z_fwf * area   ???  it is right?  I think so....
175            z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall )
176            !                                                  ! weight to respect erp field 2D structure
177            zsum_erp   = glob_sum( 'sbcfwb', ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
178            z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
179            !                                                  ! final correction term to apply
180            zerp_cor(:,:) = -1. * z_fwf_nsrf * zsurf_tospread * z_wgt(:,:)
181            !
182!!gm   ===>>>>  lbc_lnk should be useless as all the computation is done over the whole domain !
183            CALL lbc_lnk( 'sbcfwb', zerp_cor, 'T', 1. )
184            !
185            emp(:,:) = emp(:,:) + zerp_cor(:,:)
186            qns(:,:) = qns(:,:) - zerp_cor(:,:) * rcp * sst_m(:,:)  ! account for change to the heat budget due to fw correction
187            erp(:,:) = erp(:,:) + zerp_cor(:,:)
188            !
189            IF( nprint == 1 .AND. lwp ) THEN                   ! control print
190               IF( z_fwf < 0._wp ) THEN
191                  WRITE(numout,*)'   z_fwf < 0'
192                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
193               ELSE
194                  WRITE(numout,*)'   z_fwf >= 0'
195                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
196               ENDIF
197               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv'
198               WRITE(numout,*)'   z_fwf         = ', z_fwf      ,' Kg/m2/s'
199               WRITE(numout,*)'   z_fwf_nsrf    = ', z_fwf_nsrf ,' Kg/m2/s'
200               WRITE(numout,*)'   MIN(zerp_cor) = ', MINVAL(zerp_cor) 
201               WRITE(numout,*)'   MAX(zerp_cor) = ', MAXVAL(zerp_cor) 
202            ENDIF
203         ENDIF
204         DEALLOCATE( ztmsk_neg , ztmsk_pos , ztmsk_tospread , z_wgt , zerp_cor )
205         !
206      CASE DEFAULT                           !==  you should never be there  ==!
207         CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1, 2 or 3' )
208         !
209      END SELECT
210      !
211   END SUBROUTINE sbc_fwb
212
213   !!======================================================================
214END MODULE sbcfwb
Note: See TracBrowser for help on using the repository browser.