New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcfwb.F90 in trunk/NEMO/OPA_SRC/SBC – NEMO

source: trunk/NEMO/OPA_SRC/SBC/sbcfwb.F90 @ 1553

Last change on this file since 1553 was 1553, checked in by ctlod, 15 years ago

add an option to control the global freshwater budget emp in sbcfwb.F90, see ticket: #501
  • Property svn:keywords set to Id
File size: 10.2 KB
Line 
1MODULE sbcfwb
2   !!======================================================================
3   !!                       ***  MODULE  sbcfwb  ***
4   !! Ocean fluxes   : domain averaged freshwater budget
5   !!======================================================================
6   !! History :  8.2  !  01-02  (E. Durand)  Original code
7   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
8   !!            9.0  !  06-08  (G. Madec)  Surface module
9   !!            9.2  !  09-07  (C. Talandier) emp mean s spread over erp area
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   sbc_fwb      : freshwater budget for global ocean configurations
14   !!                  in free surface and 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 phycst          ! physical constants
20   USE sbcrnf          ! ocean runoffs
21   USE sbcssr          ! SS damping terms
22   USE daymod          ! calendar
23   USE in_out_manager  ! I/O manager
24   USE lib_mpp         ! distribued memory computing library
25   USE lbclnk          ! ocean lateral boundary conditions
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   sbc_fwb      ! routine called by step
31
32   REAL(wp) ::   a_fwb_b            ! annual domain averaged freshwater budget
33   REAL(wp) ::   a_fwb              ! for 2 year before (_b) and before year.
34   REAL(wp) ::   empold             ! empold to be suppressed
35   REAL(wp) ::   area               ! global mean ocean surface (interior domain)
36
37   REAL(wp), DIMENSION(jpi,jpj) ::   e1e2_i    ! area of the interior domain (e1t*e2t*tmask_i)
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !!  OPA 9.0 , LOCEAN-IPSL (2006)
44   !! $Id$
45   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE sbc_fwb( kt, kn_fwb, kn_fsbc )
50      !!---------------------------------------------------------------------
51      !!                  ***  ROUTINE sbc_fwb  ***
52      !!
53      !! ** Purpose :   Control the mean sea surface drift
54      !!
55      !! ** Method  :   several ways  depending on kn_fwb
56      !!                =0 no control
57      !!                =1 global mean of emp set to zero at each nn_fsbc time step
58      !!                =2 annual global mean corrected from previous year
59      !!                =3 global mean of emp set to zero at each nn_fsbc time step
60      !!                   & spread out over erp area depending its sign
61      !!----------------------------------------------------------------------
62      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
63      INTEGER, INTENT( in ) ::   kn_fsbc  !
64      INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index
65      !!
66      INTEGER  ::   inum                  ! temporary logical unit
67      INTEGER  ::   ikty, iyear           !
68      CHARACTER (len=32) ::   clname 
69      REAL(wp) ::   z_emp, z_emp_nsrf       ! temporary scalars
70      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread
71      REAL(wp), DIMENSION(jpi,jpj) ::   ztmsk_neg, ztmsk_pos, ztmsk_tospread
72      REAL(wp), DIMENSION(jpi,jpj) ::   z_wgt, zerp_cor
73      !!----------------------------------------------------------------------
74      !
75      IF( kt == nit000 ) THEN
76         !
77         IF(lwp) THEN
78            WRITE(numout,*)
79            WRITE(numout,*) 'sbc_fwb : FreshWater Budget correction'
80            WRITE(numout,*) '~~~~~~~'
81            IF( kn_fwb == 1 )   WRITE(numout,*) '          instantaneously set to zero'
82            IF( kn_fwb == 2 )   WRITE(numout,*) '          adjusted from previous year budget'
83            IF( kn_fwb == 3 )   WRITE(numout,*) '          set to zero and spread out over erp area'
84         ENDIF
85         !
86         e1e2_i(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
87         area = SUM( e1e2_i(:,:) )
88         IF( lk_mpp )   CALL  mpp_sum( area    )   ! sum over the global domain
89         !
90      ENDIF
91     
92
93      SELECT CASE ( kn_fwb )
94      !
95      CASE ( 0 )                               
96         WRITE(ctmp1,*)'sbc_fwb : nn_fwb = ', kn_fwb, ' is not yet associated to an option, choose either 1/2'
97         CALL ctl_stop( ctmp1 )
98         !
99         
100      !
101      CASE ( 1 )                               ! global mean emp set to zero
102         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
103            z_emp = SUM( e1e2_i(:,:) * emp(:,:) ) / area
104            IF( lk_mpp )   CALL  mpp_sum( z_emp    )   ! sum over the global domain
105            emp (:,:) = emp (:,:) - z_emp
106            emps(:,:) = emps(:,:) - z_emp
107         ENDIF
108         !
109      CASE ( 2 )                               ! emp budget adjusted from the previous year
110         ! initialisation
111         IF( kt == nit000 ) THEN
112            ! Read the corrective factor on precipitations (empold)
113            clname = 'EMPave_old.dat'
114            CALL ctlopn( inum, clname, 'OLD', 'FORMATTED', 'SEQUENTIAL',   &
115               &           1, numout, .FALSE., 1 )
116
117            READ ( inum, "(24X,I8,2ES24.16)" ) iyear, a_fwb_b, a_fwb
118            CLOSE( inum )
119            empold = a_fwb                  ! current year freshwater budget correction
120            !                               ! estimate from the previous year budget
121            IF(lwp)WRITE(numout,*)
122            IF(lwp)WRITE(numout,*)'sbc_fwb : year = ',iyear  , ' freshwater budget correction = ', empold
123            IF(lwp)WRITE(numout,*)'          year = ',iyear-1, ' freshwater budget read       = ', a_fwb
124            IF(lwp)WRITE(numout,*)'          year = ',iyear-2, ' freshwater budget read       = ', a_fwb_b
125         ENDIF   
126         !
127         ! Update empold if new year start
128         ikty = 365 * 86400 / rdttra(1)    !!bug  use of 365 days leap year or 360d year !!!!!!!
129         IF( MOD( kt, ikty ) == 0 ) THEN
130            a_fwb_b = a_fwb
131            a_fwb   = SUM( e1e2_i(:,:) * sshn(:,:) )
132            IF( lk_mpp )   CALL  mpp_sum( a_fwb    )   ! sum over the global domain
133            a_fwb   = a_fwb * 1.e+3 / ( area * 86400. * 365. )     ! convert in Kg/m3/s = mm/s
134!!gm        !                                                      !!bug 365d year
135            empold =  a_fwb                 ! current year freshwater budget correction
136            !                               ! estimate from the previous year budget
137         ENDIF
138         !
139         ! correct the freshwater fluxes
140         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
141            emp (:,:) = emp (:,:) + empold
142            emps(:,:) = emps(:,:) + empold
143         ENDIF
144         !
145         ! save empold value in a file
146         IF( kt == nitend .AND. lwp ) THEN
147            clname = 'EMPav.dat'   
148            CALL ctlopn( inum, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
149               &           1, numout, .FALSE., 0 )
150            WRITE(inum, "(24X,I8,2ES24.16)" ) nyear, a_fwb_b, a_fwb
151         ENDIF
152         !
153      CASE ( 3 )                               ! global emp set to zero and spread out over erp area
154         !
155         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
156            ! Select <0 and >0 area of erp
157            ztmsk_pos(:,:) = tmask_i(:,:)
158            WHERE( erp < 0.e0 ) ztmsk_pos = 0.e0
159            ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
160
161            ! Area filled by <0 and >0 erp
162            zsurf_neg = SUM( e1e2_i(:,:)*ztmsk_neg(:,:) )
163            zsurf_pos = SUM( e1e2_i(:,:)*ztmsk_pos(:,:) )
164       
165            ! emp global mean
166            z_emp = SUM( e1e2_i(:,:) * emp(:,:) ) / area
167            !
168            IF( lk_mpp )   CALL  mpp_sum( z_emp )
169           
170            IF( z_emp < 0.e0 ) THEN
171                ! to spread out over >0 erp area to increase evaporation damping process
172                zsurf_tospread = zsurf_pos
173                ztmsk_tospread(:,:) = ztmsk_pos(:,:)
174            ELSE
175                ! to spread out over <0 erp area to increase precipitation damping process
176                zsurf_tospread = zsurf_neg
177                ztmsk_tospread(:,:) = ztmsk_neg(:,:)
178            ENDIF
179
180            ! emp global mean over <0 or >0 erp area
181            z_emp_nsrf = SUM( e1e2_i(:,:) * z_emp ) / ( zsurf_tospread + rsmall )
182            ! weight to respect erp field 2D structure
183            z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( SUM( ztmsk_tospread(:,:) * erp(:,:) * e1e2_i(:,:) ) + rsmall )
184            ! final correction term to apply
185            zerp_cor(:,:) = -1. * z_emp_nsrf * zsurf_tospread * z_wgt(:,:)
186
187            CALL lbc_lnk( zerp_cor, 'T', 1. )
188
189            emp (:,:) = emp (:,:) + zerp_cor(:,:)
190            emps(:,:) = emps(:,:) + zerp_cor(:,:)
191            erp (:,:) = erp (:,:) + zerp_cor(:,:)
192           
193            IF( nprint == 1 .AND. lwp ) THEN
194               IF( z_emp < 0.e0 ) THEN
195                  WRITE(numout,*)'       z_emp < 0'
196                  WRITE(numout,*)'       SUM(erp+)        = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2_i(:,:) )*1.e-3,' m3.s-1'
197               ELSE
198                   WRITE(numout,*)'      z_emp >= 0'
199                   WRITE(numout,*)'      SUM(erp-)        = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2_i(:,:) )*1.e-3,' m3.s-1'
200               ENDIF
201               WRITE(numout,*)'      SUM(empG)        = ', SUM( z_emp*e1e2_i(:,:) )*1.e-3,' m3.s-1'
202               WRITE(numout,*)'      z_emp            = ', z_emp      ,' mm.s-1'
203               WRITE(numout,*)'      z_emp_nsrf       = ', z_emp_nsrf ,' mm.s-1'
204               WRITE(numout,*)'      MIN(zerp_cor)    = ', MINVAL(zerp_cor) 
205               WRITE(numout,*)'      MAX(zerp_cor)    = ', MAXVAL(zerp_cor) 
206            ENDIF
207            !
208         ENDIF
209         !
210      CASE DEFAULT    ! you should never be there
211         WRITE(ctmp1,*)'sbc_fwb : nn_fwb = ', kn_fwb, ' is not permitted for the FreshWater Budget correction, choose either 1/2'
212         CALL ctl_stop( ctmp1 )
213         !
214      END SELECT
215      !
216   END SUBROUTINE sbc_fwb
217
218   !!======================================================================
219END MODULE sbcfwb
Note: See TracBrowser for help on using the repository browser.