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.
flxfwb.F90 in trunk/NEMO/OPA_SRC/SBC – NEMO

source: trunk/NEMO/OPA_SRC/SBC/flxfwb.F90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.6 KB
Line 
1MODULE flxfwb
2   !!======================================================================
3   !!                       ***  MODULE  flxfwb  ***
4   !! Ocean fluxes   : domain averaged freshwater budget
5   !!======================================================================
6#if ! defined key_dynspg_rl
7   !!----------------------------------------------------------------------
8   !!   NOT 'key_dynspg_rl'                        Free surface formulation
9   !!----------------------------------------------------------------------
10   !!   flx_fwb      : freshwater budget for global ocean configurations
11   !!                  in free surface and forced mode
12   !!   flx_fwb_init : freshwater budget for global ocean configurations
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE cpl_oce         ! coupled atmosphere/ocean
18   USE phycst          ! physical constants
19   USE in_out_manager  ! I/O manager
20   USE lib_mpp         ! distribued memory computing library
21   USE flxrnf          ! ocean runoffs
22   USE ocesbc          ! ocean surface boudaries conditions
23   USE blk_oce
24   USE flxblk          ! bulk formulea
25   USE daymod          ! calendar
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Routine accessibility
31   PUBLIC flx_fwb      ! routine called by step
32   PUBLIC flx_fwb_init ! routine called by opa
33
34   !! * Share module variables
35   LOGICAL, PUBLIC ::    & !!: * namelist *
36      ln_fwb = .TRUE.       !: Flag to activate the fwb computation
37   REAL(wp), PUBLIC ::   &  !:
38      a_fwb_b  ,         &  !: annual domain averaged freshwater budget
39      a_fwb                 !: for 2 year before (_b) and before year.
40
41   !! * Module variables
42   REAL(wp) ::   &
43      a_emp   ,  & ! domain averaged evaporation minus precipitation
44      a_precip,  & ! domain averaged precipitation
45      a_rnf   ,  & ! domain averaged runoff
46      a_sshb  ,  & ! domain averaged sea surface heigh at nit000
47      a_sshend,  & ! domain averaged sea surface heigh at nitend
48      a_sal000,  & ! domain averaged ocean salinity at nit000 (before time step)
49      a_salend,  & ! domain averaged ocean salinity at nitend (now time step)
50      a_aminus,  & !
51      a_aplus
52      REAL(wp), DIMENSION(jpi,jpj) ::  &
53         e1e2_i    ! area of the interior domain (e1t*e2t*tmask_i)
54
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58   !!----------------------------------------------------------------------
59   !!----------------------------------------------------------------------
60   !!  OPA 9.0 , LOCEAN-IPSL (2005)
61   !! $Header$
62   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE flx_fwb( kt )
67      !!---------------------------------------------------------------------
68      !!                  ***  ROUTINE flx_fwb  ***
69      !!
70      !! ** Purpose :
71      !!
72      !! ** Method  :
73      !!
74      !! History :
75      !!   8.2  !  01-02  (E. Durand)  Original code
76      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
77      !!----------------------------------------------------------------------
78      !! * Arguments
79      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
80
81      !! * Local declarations
82      CHARACTER (len=32) ::   &
83         clname = 'EMPav.dat'
84      INTEGER  ::   ji, jj, jk        ! dummy loop indices
85      INTEGER  ::   inum              ! temporary logical unit
86      INTEGER  ::   ikty              !
87      REAL(wp) ::   &
88         zarea, zvol, zwei,       &  ! temporary scalars
89         zsm0, zempnew                !    "         "
90      !!----------------------------------------------------------------------
91
92      ! Mean global salinity
93      zsm0 = 34.72654
94
95      ! To compute emp mean value mean emp
96
97      IF( kt == nit000 ) THEN
98         IF(lwp) THEN
99            WRITE(numout,*)
100            WRITE(numout,*) 'flx_fwb : FreshWater Budget correction'
101            WRITE(numout,*) '~~~~~~~'
102         ENDIF
103
104         a_emp    = 0.e0
105         a_precip = 0.e0
106         a_rnf    = 0.e0
107         a_sshb   = 0.e0   ! averaged sea surface heigh at nit000
108         a_sal000 = 0.e0   ! averaged ocean salinity at nit000
109         a_aminus = 0.e0
110         a_aplus  = 0.e0
111         
112         e1e2_i(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
113         a_sshb = SUM( e1e2_i(:,:) * sshn(:,:) )
114         IF( lk_mpp )   CALL  mpp_sum( a_sshb   )   ! sum over all the global domain
115         DO jk = 1, jpkm1
116            DO jj = 2, jpjm1
117               DO ji = fs_2, fs_jpim1   ! vector opt.
118                  zwei  = e1e2_i(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)
119                  a_sal000 = a_sal000 + ( sb(ji,jj,jk) - zsm0 ) * zwei
120               END DO
121            END DO
122         END DO
123         IF( lk_mpp )   CALL  mpp_sum( a_sal000 )   ! sum over the global domain
124
125      ENDIF
126     
127      ! cumulate surface freshwater budget at each time-step
128      ! --------------------------------------====----------
129      a_emp    = SUM( e1e2_i(:,:) * emp   (:,:) )
130      IF( lk_mpp )   CALL  mpp_sum( a_emp    )   ! sum over the global domain
131#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
132      a_precip = SUM( e1e2_i(:,:) * watm  (:,:) )
133      IF( lk_mpp )   CALL  mpp_sum( a_precip )   ! sum over the global domain
134#endif
135      a_rnf    = SUM( e1e2_i(:,:) * runoff(:,:) )
136      IF( lk_mpp )   CALL  mpp_sum( a_rnf    )   ! sum over the global domain
137
138      IF( aminus /= 0.e0 ) a_aminus = a_aminus + ( MIN( aplus, aminus ) / aminus )
139      IF( aplus  /= 0.e0 ) a_aplus  = a_aplus  + ( MIN( aplus, aminus ) / aplus  )
140
141
142      ! Update empold if new year start
143      ikty = 365 * 86400 / rdttra(1)    !!bug  use of 365 days leap year or 360d year !!!!!!!
144      IF( MOD( kt, ikty ) == 0 ) THEN
145         zarea    = SUM( e1e2_i(:,:)            )
146         IF( lk_mpp )   CALL  mpp_sum( zarea    )   ! sum over the global domain
147         a_fwb_b = a_fwb
148         a_fwb   = SUM( e1e2_i(:,:) * sshn(:,:) )
149         IF( lk_mpp )   CALL  mpp_sum( a_fwb    )   ! sum over the global domain
150
151         a_fwb   = a_fwb * 1.e+3 / ( zarea * 86400. * 365. )    ! convert in Kg/m3/s = mm/s
152         !                                                      !!bug 365d year
153         empold =  a_fwb                 ! current year freshwater budget correction
154         !                               ! estimate from the previous year budget
155      ENDIF
156
157
158      IF( kt == nitend ) THEN
159         zvol  = 0.e0
160         zempnew = 0.e0
161         ! Mean sea level at nitend
162         a_sshend = SUM( e1e2_i(:,:) * sshn(:,:) )
163         zarea    = SUM( e1e2_i(:,:)             )
164         
165         a_salend = 0.e0
166         DO jk = 1, jpkm1   
167            DO jj = 2, jpjm1
168               DO ji = fs_2, fs_jpim1   ! vector opt.
169                  zwei  = e1e2_i(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)
170                  a_salend = a_salend + ( sn(ji,jj,jk) - zsm0 ) * zwei
171                  zvol  = zvol  + zwei
172               END DO
173            END DO
174         END DO
175         IF( lk_mpp ) THEN
176            CALL  mpp_sum( zarea    )   ! sums over all the global domain
177            CALL  mpp_sum( zvol     )       
178            CALL  mpp_sum( a_sshend )       
179            CALL  mpp_sum( a_salend )       
180         ENDIF
181
182         a_aminus = a_aminus / ( nitend - nit000 + 1 )
183         a_aplus  = a_aplus  / ( nitend - nit000 + 1 )
184
185         ! Conversion in m3
186         a_emp    = a_emp    * rdttra(1) * 1.e-3 
187         a_precip = a_precip * rdttra(1) * 1.e-3 / rday
188         a_rnf    = a_rnf    * rdttra(1) * 1.e-3
189         
190      ENDIF
191
192
193      ! Ecriture des diagnostiques
194      ! --------------------------
195
196      IF( kt == nitend ) THEN
197
198         CALL ctlopn( inum, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
199            &           1, numout, .FALSE., 0 )
200         WRITE(inum, "(24X,I8,2ES24.16)" ) nyear, a_fwb_b, a_fwb
201         WRITE(inum,*)
202         WRITE(inum,*)    'Net freshwater budget '
203         WRITE(inum,9010) '  emp    = ', a_emp   , ' m3 =', a_emp   /((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
204         WRITE(inum,9010) '  precip = ', a_precip, ' m3 =', a_precip/((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
205         WRITE(inum,9010) '  a_rnf  = ', a_rnf   , ' m3 =', a_rnf   /((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
206         WRITE(inum,*)
207         WRITE(inum,9010) '  zarea =',zarea
208         WRITE(inum,9010) '  zvol  =',zvol
209         WRITE(inum,*)
210         WRITE(inum,*)    'Mean sea level : '
211         WRITE(inum,9010) '  at nit000 = ',a_sshb           ,' m3 '
212         WRITE(inum,9010) '  at nitend = ',a_sshend         ,' m3 '
213         WRITE(inum,9010) '  diff      = ',(a_sshend-a_sshb),' m3 =', (a_sshend-a_sshb)/((nitend-nit000+1)*rdt) * 1.e-6,' Sv'
214         WRITE(inum,9020) '  mean sea level elevation    =', a_sshend/zarea,' m'
215         WRITE(inum,*)
216         WRITE(inum,*)    'Anomaly of salinity content : '
217         WRITE(inum,9010) '  at nit000 = ', a_sal000           ,' psu.m3 '
218         WRITE(inum,9010) '  at nitend = ', a_salend           ,' psu.m3 '
219         WRITE(inum,9010) '  diff      = ', a_salend - a_sal000,' psu.m3'
220         WRITE(inum,*)
221         WRITE(inum,*)    'Mean salinity : '
222         WRITE(inum,9020) '  at nit000 = ',  a_sal000/zvol+zsm0     ,' psu '
223         WRITE(inum,9020) '  at nitend = ',  a_salend/zvol+zsm0     ,' psu '
224         WRITE(inum,9020) '  diff      = ', (a_salend-a_sal000)/zvol,' psu'
225         WRITE(inum,9020) '  S-SLevitus= ',  a_salend          /zvol,' psu'
226         WRITE(inum,*)
227         WRITE(inum,*)    'Coeff : '
228         WRITE(inum,9030) '  Alpha+   =  ', a_aplus
229         WRITE(inum,9030) '  Alpha-   =  ', a_aminus
230         WRITE(inum,*)
231      ENDIF
232
233 9006 FORMAT(1X,A,ES24.16)
234 9010 FORMAT(1X,A,ES12.5,A,F10.5,A)
235 9020 FORMAT(1X,A,F10.5,A)
236 9030 FORMAT(1X,A,F8.2,A)
237
238   END SUBROUTINE flx_fwb
239
240
241   SUBROUTINE flx_fwb_init
242      !!---------------------------------------------------------------------
243      !!                ***  ROUTINE flx_fwb_init  ***
244      !!
245      !! ** Purpose :
246      !!
247      !! ** Method  :
248      !! 
249      !! History :
250      !!   9.0  !  03-09  (G. Madec)  Original code
251      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
252      !!----------------------------------------------------------------------
253      !! * Local declarations
254      CHARACTER (len=32) ::   &
255         clname = 'EMPave_old.dat'
256      INTEGER ::   inum, iyear
257
258      NAMELIST/namfwb/ ln_fwb 
259      !!----------------------------------------------------------------------
260         
261      ! Read Namelist namfwb : freshWater Budget correction
262      ! --------------------
263      REWIND( numnam )
264      READ  ( numnam, namfwb )           
265     
266      ! Parameter control and print
267      ! ---------------------------
268      ! Control print
269      IF(lwp) THEN
270         WRITE(numout,*)
271         WRITE(numout,*) 'flx_fwb_init : FreshWater Budget correction'
272         WRITE(numout,*) '~~~~~~~~~~~~'
273         WRITE(numout,*) '               Namelist namfwb : set fwb parameters'
274         WRITE(numout,*) '                  use or not fwb correction      ln_fwb   = ', ln_fwb
275      ENDIF
276      ! Option consistency
277#if defined key_dynspg_rl
278      IF(lwp) WRITE '               Rigid-lid option, fwb correction is useless, but valid'
279#endif
280      IF( lk_cpl ) THEN
281         IF(lwp) WRITE(numout,*) '               Coupled option, fwb correction is a flux correction ! '
282         IF(lwp) WRITE(numout,*) '               ln_fwb = .FALSE. is recommanded'
283      ENDIF
284
285      IF( nbit_cmp == 1 .AND. ln_fwb ) THEN
286         CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require ln_fwb = .false.' )
287      END IF
288
289      !                                        ! ==============================
290      IF( ln_fwb ) THEN                        !  Freshwater budget correction
291         !                                     ! ==============================
292         ! Read the corrective factor on precipitations (empold)
293         CALL ctlopn( inum, clname, 'OLD', 'FORMATTED', 'SEQUENTIAL',   &
294            &           1, numout, .FALSE., 1 )
295         READ ( inum, "(24X,I8,2ES24.16)" ) iyear, a_fwb_b, a_fwb
296         CLOSE( inum )
297         empold = a_fwb                  ! current year freshwater budget correction
298         !                               ! estimate from the previous year budget
299         IF(lwp)WRITE(numout,*)
300         IF(lwp)WRITE(numout,*)'flx_fwb_init : year = ',iyear  , ' freshwater budget correction = ', empold
301         IF(lwp)WRITE(numout,*)'               year = ',iyear-1, ' freshwater budget read       = ', a_fwb
302         IF(lwp)WRITE(numout,*)'               year = ',iyear-2, ' freshwater budget read       = ', a_fwb_b
303         !                                    ! ==============================
304      ELSE                                    !      NO  budget correction
305         !                                    ! ==============================
306         IF(lwp) WRITE(numout,*)
307         IF(lwp) WRITE(numout,*)'               NO freshwater budget correction'
308         empold = 0.e0
309      ENDIF
310
311   END SUBROUTINE flx_fwb_init
312
313#else
314   !!----------------------------------------------------------------------
315   !!   Default case :                       
316   !!----------------------------------------------------------------------
317   LOGICAL, PUBLIC ::   ln_fwb = .FALSE.   !: no fwb forced
318CONTAINS
319   SUBROUTINE flx_fwb( kt )                ! dummy routine
320      WRITE(*,*) 'flx_fwb: You should not have seen this print! error?', kt
321   END SUBROUTINE flx_fwb
322   SUBROUTINE flx_fwb_init
323   END SUBROUTINE flx_fwb_init
324#endif
325   !!======================================================================
326END MODULE flxfwb
Note: See TracBrowser for help on using the repository browser.