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
RevLine 
[888]1MODULE sbcfwb
2   !!======================================================================
3   !!                       ***  MODULE  sbcfwb  ***
4   !! Ocean fluxes   : domain averaged freshwater budget
5   !!======================================================================
[2528]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
[888]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
[1553]21   USE sbcssr          ! SS damping terms
[888]22   USE in_out_manager  ! I/O manager
23   USE lib_mpp         ! distribued memory computing library
[3294]24   USE wrk_nemo        ! work arrays
25   USE timing          ! Timing
[1553]26   USE lbclnk          ! ocean lateral boundary conditions
[2528]27   USE lib_fortran
[888]28
29   IMPLICIT NONE
30   PRIVATE
31
[2715]32   PUBLIC   sbc_fwb    ! routine called by step
[888]33
[2715]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)
[3964]38   REAL(wp) ::   r1_rau0   ! = 1 / rau0
[888]39
40   !! * Substitutions
41#  include "domzgr_substitute.h90"
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
[2528]44   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1156]45   !! $Id$
[2528]46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[888]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
[1168]58      !!                =1 global mean of emp set to zero at each nn_fsbc time step
59      !!                =2 annual global mean corrected from previous year
[1553]60      !!                =3 global mean of emp set to zero at each nn_fsbc time step
61      !!                   & spread out over erp area depending its sign
[3962]62      !! Note: if mass exchanges between ice and ocean, it is taken into account
63      !!       when computing the budget
[888]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
[2715]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    !   -      -
[3294]72      REAL(wp), POINTER, DIMENSION(:,:) ::   ztmsk_neg, ztmsk_pos, ztmsk_tospread, z_wgt, zerp_cor
[888]73      !!----------------------------------------------------------------------
74      !
[3294]75      IF( nn_timing == 1 )  CALL timing_start('sbc_fwb')
[2715]76      !
[3294]77      CALL wrk_alloc( jpi,jpj, ztmsk_neg, ztmsk_pos, ztmsk_tospread, z_wgt, zerp_cor )
78      !
[888]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'
[2528]86            IF( kn_fwb == 3 )   WRITE(numout,*) '          fwf set to zero and spread out over erp area'
[888]87         ENDIF
88         !
[2471]89         IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' )
90         !
[2715]91         area = glob_sum( e1e2t(:,:) )           ! interior global domain surface
[3962]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         !
[888]100      ENDIF
101     
102
103      SELECT CASE ( kn_fwb )
104      !
[2528]105      CASE ( 1 )                             !==  global mean fwf set to zero  ==!
[1245]106         !
[888]107         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
[3962]108            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - snwice_fmass(:,:) ) ) / area   ! sum over the global domain
[2528]109            emp (:,:) = emp (:,:) - z_fwf 
[3938]110            erp (:,:) = erp (:,:) - z_fwf 
[3952]111#if defined key_lim2
112            emps(:,:) = emps(:,:) - z_fwf      ! emps is the mass flux in lim2 case
113#endif
[888]114         ENDIF
115         !
[2528]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)
[1581]120            CALL ctl_opn( inum, 'EMPave_old.dat', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
[895]121            READ ( inum, "(24X,I8,2ES24.16)" ) iyear, a_fwb_b, a_fwb
122            CLOSE( inum )
[2528]123            fwfold = a_fwb                            ! current year freshwater budget correction
124            !                                         ! estimate from the previous year budget
[895]125            IF(lwp)WRITE(numout,*)
[2528]126            IF(lwp)WRITE(numout,*)'sbc_fwb : year = ',iyear  , ' freshwater budget correction = ', fwfold
[895]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   
[2528]130         !                                         ! Update fwfold if new year start
[888]131         ikty = 365 * 86400 / rdttra(1)    !!bug  use of 365 days leap year or 360d year !!!!!!!
132         IF( MOD( kt, ikty ) == 0 ) THEN
[3964]133            !
134            r1_rau0 = 1._wp / rau0
135            !
[3962]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 ) )
[888]139            a_fwb   = a_fwb * 1.e+3 / ( area * 86400. * 365. )     ! convert in Kg/m3/s = mm/s
140!!gm        !                                                      !!bug 365d year
[2528]141            fwfold =  a_fwb                           ! current year freshwater budget correction
142            !                                         ! estimate from the previous year budget
[888]143         ENDIF
[3977]144         IF( MOD( kt, 5475 ) == 0 ) THEN
145            IF(lwp) WRITE(numout,*) 'correction of fresh water budget :', fwfold
146         ENDIF 
[888]147         !
[2528]148         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN      ! correct the freshwater fluxes
149            emp (:,:) = emp (:,:) + fwfold
[3938]150            erp (:,:) = erp (:,:) + fwfold
[3952]151#if defined key_lim2
152            emps(:,:) = emps(:,:) + fwfold      ! emps is the mass flux in lim2 case
153#endif
[888]154         ENDIF
155         !
[2528]156         IF( kt == nitend .AND. lwp ) THEN         ! save fwfold value in a file
[1581]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 )
[888]160         ENDIF
161         !
[2528]162      CASE ( 3 )                             !==  global fwf set to zero and spread out over erp area  ==!
[1553]163         !
164         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
[2528]165            ztmsk_pos(:,:) = tmask_i(:,:)                      ! Select <0 and >0 area of erp
166            WHERE( erp < 0._wp )   ztmsk_pos = 0._wp
[1553]167            ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
168            !
[2715]169            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )   ! Area filled by <0 and >0 erp
170            zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
[3962]171            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges)
172            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - snwice_fmass(:,:) ) ) / area
[2528]173            !           
174            IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation
175                zsurf_tospread      = zsurf_pos
[1553]176                ztmsk_tospread(:,:) = ztmsk_pos(:,:)
[2528]177            ELSE                             ! spread out over <0 erp area to increase precipitation
178                zsurf_tospread      = zsurf_neg
[1553]179                ztmsk_tospread(:,:) = ztmsk_neg(:,:)
180            ENDIF
[2528]181            !
[2715]182            zsum_fwf   = glob_sum( e1e2t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
[2528]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
[3962]186            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
[1822]187            z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
[2528]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 !
[1553]192            CALL lbc_lnk( zerp_cor, 'T', 1. )
[2528]193            !
[1553]194            emp (:,:) = emp (:,:) + zerp_cor(:,:)
[3952]195#if defined key_lim2
196            emps(:,:) = emps(:,:) + zerp_cor(:,:)      ! emps is the mass flux in lim2 case
197#endif
[1553]198            erp (:,:) = erp (:,:) + zerp_cor(:,:)
[2528]199            !
200            IF( nprint == 1 .AND. lwp ) THEN                   ! control print
201               IF( z_fwf < 0._wp ) THEN
202                  WRITE(numout,*)'   z_fwf < 0'
[2715]203                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
[1553]204               ELSE
[2528]205                  WRITE(numout,*)'   z_fwf >= 0'
[2715]206                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
[1553]207               ENDIF
[2715]208               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv'
[2528]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) 
[1553]213            ENDIF
214         ENDIF
215         !
[2528]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' )
[888]218         !
219      END SELECT
220      !
[3294]221      CALL wrk_dealloc( jpi,jpj, ztmsk_neg, ztmsk_pos, ztmsk_tospread, z_wgt, zerp_cor )
[2715]222      !
[3294]223      IF( nn_timing == 1 )  CALL timing_stop('sbc_fwb')
224      !
[888]225   END SUBROUTINE sbc_fwb
226
227   !!======================================================================
228END MODULE sbcfwb
Note: See TracBrowser for help on using the repository browser.