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 branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90 @ 3977

Last change on this file since 3977 was 3977, checked in by flavoni, 11 years ago

add print of fwb correction value, and hminrhg for LIM2, in CNRS LIM3 branch

  • Property svn:keywords set to Id
File size: 11.6 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   !!----------------------------------------------------------------------
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 in_out_manager  ! I/O manager
23   USE lib_mpp         ! distribued memory computing library
24   USE wrk_nemo        ! work arrays
25   USE timing          ! Timing
26   USE lbclnk          ! ocean lateral boundary conditions
27   USE lib_fortran
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   sbc_fwb    ! routine called by step
33
34   REAL(wp) ::   a_fwb_b   ! annual domain averaged freshwater budget
35   REAL(wp) ::   a_fwb     ! for 2 year before (_b) and before year.
36   REAL(wp) ::   fwfold    ! fwfold to be suppressed
37   REAL(wp) ::   area      ! global mean ocean surface (interior domain)
38   REAL(wp) ::   r1_rau0   ! = 1 / rau0
39
40   !! * Substitutions
41#  include "domzgr_substitute.h90"
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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 mass exchanges between ice and ocean, it is taken into account
63      !!       when computing the budget
64      !!----------------------------------------------------------------------
65      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
66      INTEGER, INTENT( in ) ::   kn_fsbc  !
67      INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index
68      !
69      INTEGER  ::   inum, ikty, iyear   ! local integers
70      REAL(wp) ::   z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp   ! local scalars
71      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread    !   -      -
72      REAL(wp), POINTER, DIMENSION(:,:) ::   ztmsk_neg, ztmsk_pos, ztmsk_tospread, z_wgt, zerp_cor
73      !!----------------------------------------------------------------------
74      !
75      IF( nn_timing == 1 )  CALL timing_start('sbc_fwb')
76      !
77      CALL wrk_alloc( jpi,jpj, ztmsk_neg, ztmsk_pos, ztmsk_tospread, z_wgt, zerp_cor )
78      !
79      IF( kt == nit000 ) THEN
80         IF(lwp) THEN
81            WRITE(numout,*)
82            WRITE(numout,*) 'sbc_fwb : FreshWater Budget correction'
83            WRITE(numout,*) '~~~~~~~'
84            IF( kn_fwb == 1 )   WRITE(numout,*) '          instantaneously set to zero'
85            IF( kn_fwb == 2 )   WRITE(numout,*) '          adjusted from previous year budget'
86            IF( kn_fwb == 3 )   WRITE(numout,*) '          fwf set to zero and spread out over erp area'
87         ENDIF
88         !
89         IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' )
90         !
91         area = glob_sum( e1e2t(:,:) )           ! interior global domain surface
92         !
93!!#if ! defined key_lim2 &&  ! defined key_lim3 && ! defined key_cice
94#if ! defined key_lim3
95         snwice_mass_b(:,:) = 0._wp              ! no sea-ice model is being used : no snow+ice mass
96         snwice_mass  (:,:) = 0._wp
97         snwice_fmass (:,:) = 0._wp 
98#endif
99         !
100      ENDIF
101     
102
103      SELECT CASE ( kn_fwb )
104      !
105      CASE ( 1 )                             !==  global mean fwf set to zero  ==!
106         !
107         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
108            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - snwice_fmass(:,:) ) ) / area   ! sum over the global domain
109            emp (:,:) = emp (:,:) - z_fwf 
110            erp (:,:) = erp (:,:) - z_fwf 
111#if defined key_lim2
112            emps(:,:) = emps(:,:) - z_fwf      ! emps is the mass flux in lim2 case
113#endif
114         ENDIF
115         !
116      CASE ( 2 )                             !==  fwf budget adjusted from the previous year  ==!
117         !
118         IF( kt == nit000 ) THEN                   ! initialisation
119            !                                         ! Read the corrective factor on precipitations (fwfold)
120            CALL ctl_opn( inum, 'EMPave_old.dat', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
121            READ ( inum, "(24X,I8,2ES24.16)" ) iyear, a_fwb_b, a_fwb
122            CLOSE( inum )
123            fwfold = a_fwb                            ! current year freshwater budget correction
124            !                                         ! estimate from the previous year budget
125            IF(lwp)WRITE(numout,*)
126            IF(lwp)WRITE(numout,*)'sbc_fwb : year = ',iyear  , ' freshwater budget correction = ', fwfold
127            IF(lwp)WRITE(numout,*)'          year = ',iyear-1, ' freshwater budget read       = ', a_fwb
128            IF(lwp)WRITE(numout,*)'          year = ',iyear-2, ' freshwater budget read       = ', a_fwb_b
129         ENDIF   
130         !                                         ! Update fwfold if new year start
131         ikty = 365 * 86400 / rdttra(1)    !!bug  use of 365 days leap year or 360d year !!!!!!!
132         IF( MOD( kt, ikty ) == 0 ) THEN
133            !
134            r1_rau0 = 1._wp / rau0
135            !
136            a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
137                                                      ! sum over the global domain
138            a_fwb   = glob_sum( e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) )
139            a_fwb   = a_fwb * 1.e+3 / ( area * 86400. * 365. )     ! convert in Kg/m3/s = mm/s
140!!gm        !                                                      !!bug 365d year
141            fwfold =  a_fwb                           ! current year freshwater budget correction
142            !                                         ! estimate from the previous year budget
143         ENDIF
144         IF( MOD( kt, 5475 ) == 0 ) THEN
145            IF(lwp) WRITE(numout,*) 'correction of fresh water budget :', fwfold
146         ENDIF 
147         !
148         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN      ! correct the freshwater fluxes
149            emp (:,:) = emp (:,:) + fwfold
150            erp (:,:) = erp (:,:) + fwfold
151#if defined key_lim2
152            emps(:,:) = emps(:,:) + fwfold      ! emps is the mass flux in lim2 case
153#endif
154         ENDIF
155         !
156         IF( kt == nitend .AND. lwp ) THEN         ! save fwfold value in a file
157            CALL ctl_opn( inum, 'EMPave.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
158            WRITE( inum, "(24X,I8,2ES24.16)" ) nyear, a_fwb_b, a_fwb
159            CLOSE( inum )
160         ENDIF
161         !
162      CASE ( 3 )                             !==  global fwf set to zero and spread out over erp area  ==!
163         !
164         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
165            ztmsk_pos(:,:) = tmask_i(:,:)                      ! Select <0 and >0 area of erp
166            WHERE( erp < 0._wp )   ztmsk_pos = 0._wp
167            ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
168            !
169            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )   ! Area filled by <0 and >0 erp
170            zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
171            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges)
172            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - snwice_fmass(:,:) ) ) / area
173            !           
174            IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation
175                zsurf_tospread      = zsurf_pos
176                ztmsk_tospread(:,:) = ztmsk_pos(:,:)
177            ELSE                             ! spread out over <0 erp area to increase precipitation
178                zsurf_tospread      = zsurf_neg
179                ztmsk_tospread(:,:) = ztmsk_neg(:,:)
180            ENDIF
181            !
182            zsum_fwf   = glob_sum( e1e2t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
183!!gm :  zsum_fwf   = z_fwf * area   ???  it is right?  I think so....
184            z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall )
185            !                                                  ! weight to respect erp field 2D structure
186            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
187            z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
188            !                                                  ! final correction term to apply
189            zerp_cor(:,:) = -1. * z_fwf_nsrf * zsurf_tospread * z_wgt(:,:)
190            !
191!!gm   ===>>>>  lbc_lnk should be useless as all the computation is done over the whole domain !
192            CALL lbc_lnk( zerp_cor, 'T', 1. )
193            !
194            emp (:,:) = emp (:,:) + zerp_cor(:,:)
195#if defined key_lim2
196            emps(:,:) = emps(:,:) + zerp_cor(:,:)      ! emps is the mass flux in lim2 case
197#endif
198            erp (:,:) = erp (:,:) + zerp_cor(:,:)
199            !
200            IF( nprint == 1 .AND. lwp ) THEN                   ! control print
201               IF( z_fwf < 0._wp ) THEN
202                  WRITE(numout,*)'   z_fwf < 0'
203                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
204               ELSE
205                  WRITE(numout,*)'   z_fwf >= 0'
206                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
207               ENDIF
208               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv'
209               WRITE(numout,*)'   z_fwf         = ', z_fwf      ,' Kg/m2/s'
210               WRITE(numout,*)'   z_fwf_nsrf    = ', z_fwf_nsrf ,' Kg/m2/s'
211               WRITE(numout,*)'   MIN(zerp_cor) = ', MINVAL(zerp_cor) 
212               WRITE(numout,*)'   MAX(zerp_cor) = ', MAXVAL(zerp_cor) 
213            ENDIF
214         ENDIF
215         !
216      CASE DEFAULT                           !==  you should never be there  ==!
217         CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1, 2 or 3' )
218         !
219      END SELECT
220      !
221      CALL wrk_dealloc( jpi,jpj, ztmsk_neg, ztmsk_pos, ztmsk_tospread, z_wgt, zerp_cor )
222      !
223      IF( nn_timing == 1 )  CALL timing_stop('sbc_fwb')
224      !
225   END SUBROUTINE sbc_fwb
226
227   !!======================================================================
228END MODULE sbcfwb
Note: See TracBrowser for help on using the repository browser.