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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90 @ 4416

Last change on this file since 4416 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 10.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 lbclnk          ! ocean lateral boundary conditions
25   USE lib_fortran
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) ::   fwfold    ! fwfold to be suppressed
35   REAL(wp) ::   area      ! global mean ocean surface (interior domain)
36
37   !! * Control permutation of array indices
38#  include "oce_ftrans.h90"
39#  include "dom_oce_ftrans.h90"
40#  include "sbc_oce_ftrans.h90"
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE sbc_fwb( kt, kn_fwb, kn_fsbc )
53      !!---------------------------------------------------------------------
54      !!                  ***  ROUTINE sbc_fwb  ***
55      !!
56      !! ** Purpose :   Control the mean sea surface drift
57      !!
58      !! ** Method  :   several ways  depending on kn_fwb
59      !!                =0 no control
60      !!                =1 global mean of emp set to zero at each nn_fsbc time step
61      !!                =2 annual global mean corrected from previous year
62      !!                =3 global mean of emp set to zero at each nn_fsbc time step
63      !!                   & spread out over erp area depending its sign
64      !!----------------------------------------------------------------------
65      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
66      USE wrk_nemo, ONLY:   ztmsk_neg      => wrk_2d_1 , ztmsk_pos => wrk_2d_2
67      USE wrk_nemo, ONLY:   ztmsk_tospread => wrk_2d_3
68      USE wrk_nemo, ONLY:   z_wgt          => wrk_2d_4 , zerp_cor  => wrk_2d_5
69      !
70      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
71      INTEGER, INTENT( in ) ::   kn_fsbc  !
72      INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index
73      !
74      INTEGER  ::   inum, ikty, iyear   ! local integers
75      REAL(wp) ::   z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp   ! local scalars
76      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread    !   -      -
77      !!----------------------------------------------------------------------
78      !
79      IF( wrk_in_use(2, 1,2,3,4,5) ) THEN
80         CALL ctl_stop('sbc_fwb: requested workspace arrays are unavailable')   ;   RETURN
81      ENDIF
82      !
83      IF( kt == nit000 ) THEN
84         IF(lwp) THEN
85            WRITE(numout,*)
86            WRITE(numout,*) 'sbc_fwb : FreshWater Budget correction'
87            WRITE(numout,*) '~~~~~~~'
88            IF( kn_fwb == 1 )   WRITE(numout,*) '          instantaneously set to zero'
89            IF( kn_fwb == 2 )   WRITE(numout,*) '          adjusted from previous year budget'
90            IF( kn_fwb == 3 )   WRITE(numout,*) '          fwf set to zero and spread out over erp area'
91         ENDIF
92         !
93         IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' )
94         !
95         area = glob_sum( e1e2t(:,:) )           ! interior global domain surface
96      ENDIF
97     
98
99      SELECT CASE ( kn_fwb )
100      !
101      CASE ( 1 )                             !==  global mean fwf set to zero  ==!
102         !
103         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
104            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area   ! sum over the global domain
105            emp (:,:) = emp (:,:) - z_fwf 
106            emps(:,:) = emps(:,:) - z_fwf 
107         ENDIF
108         !
109      CASE ( 2 )                             !==  fwf budget adjusted from the previous year  ==!
110         !
111         IF( kt == nit000 ) THEN                   ! initialisation
112            !                                         ! Read the corrective factor on precipitations (fwfold)
113            CALL ctl_opn( inum, 'EMPave_old.dat', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
114            READ ( inum, "(24X,I8,2ES24.16)" ) iyear, a_fwb_b, a_fwb
115            CLOSE( inum )
116            fwfold = a_fwb                            ! current year freshwater budget correction
117            !                                         ! estimate from the previous year budget
118            IF(lwp)WRITE(numout,*)
119            IF(lwp)WRITE(numout,*)'sbc_fwb : year = ',iyear  , ' freshwater budget correction = ', fwfold
120            IF(lwp)WRITE(numout,*)'          year = ',iyear-1, ' freshwater budget read       = ', a_fwb
121            IF(lwp)WRITE(numout,*)'          year = ',iyear-2, ' freshwater budget read       = ', a_fwb_b
122         ENDIF   
123         !                                         ! Update fwfold if new year start
124         ikty = 365 * 86400 / rdttra(1)    !!bug  use of 365 days leap year or 360d year !!!!!!!
125         IF( MOD( kt, ikty ) == 0 ) THEN
126            a_fwb_b = a_fwb
127            a_fwb   = glob_sum( e1e2t(:,:) * sshn(:,:) )   ! sum over the global domain
128            a_fwb   = a_fwb * 1.e+3 / ( area * 86400. * 365. )     ! convert in Kg/m3/s = mm/s
129!!gm        !                                                      !!bug 365d year
130            fwfold =  a_fwb                           ! current year freshwater budget correction
131            !                                         ! estimate from the previous year budget
132         ENDIF
133         !
134         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN      ! correct the freshwater fluxes
135            emp (:,:) = emp (:,:) + fwfold
136            emps(:,:) = emps(:,:) + fwfold
137         ENDIF
138         !
139         IF( kt == nitend .AND. lwp ) THEN         ! save fwfold value in a file
140            CALL ctl_opn( inum, 'EMPave.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
141            WRITE( inum, "(24X,I8,2ES24.16)" ) nyear, a_fwb_b, a_fwb
142            CLOSE( inum )
143         ENDIF
144         !
145      CASE ( 3 )                             !==  global fwf set to zero and spread out over erp area  ==!
146         !
147         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
148            ztmsk_pos(:,:) = tmask_i(:,:)                      ! Select <0 and >0 area of erp
149            WHERE( erp < 0._wp )   ztmsk_pos = 0._wp
150            ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
151            !
152            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )   ! Area filled by <0 and >0 erp
153            zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
154            !                                                  ! fwf global mean
155            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area
156            !           
157            IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation
158                zsurf_tospread      = zsurf_pos
159                ztmsk_tospread(:,:) = ztmsk_pos(:,:)
160            ELSE                             ! spread out over <0 erp area to increase precipitation
161                zsurf_tospread      = zsurf_neg
162                ztmsk_tospread(:,:) = ztmsk_neg(:,:)
163            ENDIF
164            !
165            zsum_fwf   = glob_sum( e1e2t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
166!!gm :  zsum_fwf   = z_fwf * area   ???  it is right?  I think so....
167            z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall )
168            !                                                  ! weight to respect erp field 2D structure
169            zsum_erp = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
170            z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
171            !                                                  ! final correction term to apply
172            zerp_cor(:,:) = -1. * z_fwf_nsrf * zsurf_tospread * z_wgt(:,:)
173            !
174!!gm   ===>>>>  lbc_lnk should be useless as all the computation is done over the whole domain !
175            CALL lbc_lnk( zerp_cor, 'T', 1. )
176            !
177            emp (:,:) = emp (:,:) + zerp_cor(:,:)
178            emps(:,:) = emps(:,:) + zerp_cor(:,:)
179            erp (:,:) = erp (:,:) + zerp_cor(:,:)
180            !
181            IF( nprint == 1 .AND. lwp ) THEN                   ! control print
182               IF( z_fwf < 0._wp ) THEN
183                  WRITE(numout,*)'   z_fwf < 0'
184                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
185               ELSE
186                  WRITE(numout,*)'   z_fwf >= 0'
187                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
188               ENDIF
189               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv'
190               WRITE(numout,*)'   z_fwf         = ', z_fwf      ,' Kg/m2/s'
191               WRITE(numout,*)'   z_fwf_nsrf    = ', z_fwf_nsrf ,' Kg/m2/s'
192               WRITE(numout,*)'   MIN(zerp_cor) = ', MINVAL(zerp_cor) 
193               WRITE(numout,*)'   MAX(zerp_cor) = ', MAXVAL(zerp_cor) 
194            ENDIF
195         ENDIF
196         !
197      CASE DEFAULT                           !==  you should never be there  ==!
198         CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1, 2 or 3' )
199         !
200      END SELECT
201      !
202      IF( wrk_not_released(2, 1,2,3,4,5) )   CALL ctl_stop('sbc_fwb: failed to release workspace arrays')
203      !
204   END SUBROUTINE sbc_fwb
205
206   !!======================================================================
207END MODULE sbcfwb
Note: See TracBrowser for help on using the repository browser.