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 @ 193

Last change on this file since 193 was 193, checked in by opalod, 19 years ago

CT : BUGFIX134 : correction of the current year freshwater budget empold using the previous year budget instead of the 2 previous years budget
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.2 KB
Line 
1MODULE flxfwb
2   !!======================================================================
3   !!                       ***  MODULE  flxfwb  ***
4   !! Ocean fluxes   : domain averaged freshwater budget
5   !!======================================================================
6#if defined key_dynspg_fsc
7   !!----------------------------------------------------------------------
8   !!   'key_dynspg_fsc'               No free surface with constant volume
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
60CONTAINS
61
62   SUBROUTINE flx_fwb( kt )
63      !!---------------------------------------------------------------------
64      !!                  ***  ROUTINE flx_fwb  ***
65      !!
66      !! ** Purpose :
67      !!
68      !! ** Method  :
69      !!
70      !! History :
71      !!   8.2  !  01-02  (E. Durand)  Original code
72      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
73      !!----------------------------------------------------------------------
74      !! * Arguments
75      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
76
77      !! * Local declarations
78      INTEGER  ::   ji, jj, jk        ! dummy loop indices
79      INTEGER  ::   inum = 11         ! temporary logical unit
80      INTEGER  ::   ikty              !
81      REAL(wp) ::   &
82         zarea, zvol, zwei,       &  ! temporary scalars
83         zsm0, zempnew                !    "         "
84      !!----------------------------------------------------------------------
85
86      ! Mean global salinity
87      zsm0 = 34.72654
88
89      ! To compute emp mean value mean emp
90
91      IF( kt == nit000 ) THEN
92         IF(lwp) THEN
93            WRITE(numout,*)
94            WRITE(numout,*) 'flx_fwb : FreshWater Budget correction'
95            WRITE(numout,*) '~~~~~~~'
96         ENDIF
97
98         a_emp    = 0.e0
99         a_precip = 0.e0
100         a_rnf    = 0.e0
101         a_sshb   = 0.e0   ! averaged sea surface heigh at nit000
102         a_sal000 = 0.e0   ! averaged ocean salinity at nit000
103         a_aminus = 0.e0
104         a_aplus  = 0.e0
105         
106         e1e2_i(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
107         a_sshb = SUM( e1e2_i(:,:) * sshn(:,:) )
108         IF( lk_mpp )   CALL  mpp_sum( a_sshb   )   ! sum over all the global domain
109         DO jk = 1, jpkm1
110            DO jj = 2, jpjm1
111               DO ji = fs_2, fs_jpim1   ! vector opt.
112                  zwei  = e1e2_i(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)
113                  a_sal000 = a_sal000 + ( sb(ji,jj,jk) - zsm0 ) * zwei
114               END DO
115            END DO
116         END DO
117         IF( lk_mpp )   CALL  mpp_sum( a_sal000 )   ! sum over the global domain
118
119      ENDIF
120     
121      ! cumulate surface freshwater budget at each time-step
122      ! --------------------------------------====----------
123      a_emp    = SUM( e1e2_i(:,:) * emp   (:,:) )
124      IF( lk_mpp )   CALL  mpp_sum( a_emp    )   ! sum over the global domain
125#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
126      a_precip = SUM( e1e2_i(:,:) * watm  (:,:) )
127      IF( lk_mpp )   CALL  mpp_sum( a_precip )   ! sum over the global domain
128#endif
129      a_rnf    = SUM( e1e2_i(:,:) * runoff(:,:) )
130      IF( lk_mpp )   CALL  mpp_sum( a_rnf    )   ! sum over the global domain
131
132      IF( aminus /= 0.e0 ) a_aminus = a_aminus + ( MIN( aplus, aminus ) / aminus )
133      IF( aplus  /= 0.e0 ) a_aplus  = a_aplus  + ( MIN( aplus, aminus ) / aplus  )
134
135
136      ! Update empold if new year start
137      ikty = 365 * 86400 / rdttra(1)    !!bug  use of 365 days leap year or 360d year !!!!!!!
138      IF( MOD( kt, ikty ) == 0 ) THEN
139         zarea    = SUM( e1e2_i(:,:)            )
140         IF( lk_mpp )   CALL  mpp_sum( zarea    )   ! sum over the global domain
141         a_fwb_b = a_fwb
142         a_fwb   = SUM( e1e2_i(:,:) * sshn(:,:) )
143         IF( lk_mpp )   CALL  mpp_sum( a_fwb    )   ! sum over the global domain
144
145         a_fwb   = a_fwb * 1.e+3 / ( zarea * 86400. * 365. )    ! convert in Kg/m3/s = mm/s
146         !                                                      !!bug 365d year
147         empold =  a_fwb                 ! current year freshwater budget correction
148         !                               ! estimate from the previous year budget
149      ENDIF
150
151
152      IF( kt == nitend ) THEN
153         zvol  = 0.e0
154         zempnew = 0.e0
155         ! Mean sea level at nitend
156         a_sshend = SUM( e1e2_i(:,:) * sshn(:,:) )
157         zarea    = SUM( e1e2_i(:,:)             )
158         
159         a_salend = 0.e0
160         DO jk = 1, jpkm1   
161            DO jj = 2, jpjm1
162               DO ji = fs_2, fs_jpim1   ! vector opt.
163                  zwei  = e1e2_i(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)
164                  a_salend = a_salend + ( sn(ji,jj,jk) - zsm0 ) * zwei
165                  zvol  = zvol  + zwei
166               END DO
167            END DO
168         END DO
169         IF( lk_mpp ) THEN
170            CALL  mpp_sum( zarea    )   ! sums over all the global domain
171            CALL  mpp_sum( a_sshend )       
172            CALL  mpp_sum( a_salend )       
173         ENDIF
174
175         a_aminus = a_aminus / ( nitend - nit000 + 1 )
176         a_aplus  = a_aplus  / ( nitend - nit000 + 1 )
177
178         ! Conversion in m3
179         a_emp    = a_emp    * rdttra(1) * 1.e-3 
180         a_precip = a_precip * rdttra(1) * 1.e-3 / rday
181         a_rnf    = a_rnf    * rdttra(1) * 1.e-3
182         
183      ENDIF
184
185
186      ! Ecriture des diagnostiques
187      ! --------------------------
188
189      IF( kt == nitend ) THEN
190
191         OPEN( inum, FILE='EMPave.dat' )
192         WRITE(inum, "(24X,I8,2ES24.16)" ) nyear, a_fwb_b, a_fwb
193         WRITE(inum,*)
194         WRITE(inum,*)    'Net freshwater budget '
195         WRITE(inum,9010) '  emp    = ', a_emp   , ' m3 =', a_emp   /((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
196         WRITE(inum,9010) '  precip = ', a_precip, ' m3 =', a_precip/((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
197         WRITE(inum,9010) '  a_rnf  = ', a_rnf   , ' m3 =', a_rnf   /((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
198         WRITE(inum,*)
199         WRITE(inum,9010) '  zarea =',zarea
200         WRITE(inum,9010) '  zvol  =',zvol
201         WRITE(inum,*)
202         WRITE(inum,*)    'Mean sea level : '
203         WRITE(inum,9010) '  at nit000 = ',a_sshb           ,' m3 '
204         WRITE(inum,9010) '  at nitend = ',a_sshend         ,' m3 '
205         WRITE(inum,9010) '  diff      = ',(a_sshend-a_sshb),' m3 =', (a_sshend-a_sshb)/((nitend-nit000+1)*rdt) * 1.e-6,' Sv'
206         WRITE(inum,9020) '  mean sea level elevation    =', a_sshend/zarea,' m'
207         WRITE(inum,*)
208         WRITE(inum,*)    'Anomaly of salinity content : '
209         WRITE(inum,9010) '  at nit000 = ', a_sal000           ,' psu.m3 '
210         WRITE(inum,9010) '  at nitend = ', a_salend           ,' psu.m3 '
211         WRITE(inum,9010) '  diff      = ', a_salend - a_sal000,' psu.m3'
212         WRITE(inum,*)
213         WRITE(inum,*)    'Mean salinity : '
214         WRITE(inum,9020) '  at nit000 = ',  a_sal000/zvol+zsm0     ,' psu '
215         WRITE(inum,9020) '  at nitend = ',  a_salend/zvol+zsm0     ,' psu '
216         WRITE(inum,9020) '  diff      = ', (a_salend-a_sal000)/zvol,' psu'
217         WRITE(inum,9020) '  S-SLevitus= ',  a_salend          /zvol,' psu'
218         WRITE(inum,*)
219         WRITE(inum,*)    'Coeff : '
220         WRITE(inum,9030) '  Alpha+   =  ', a_aplus
221         WRITE(inum,9030) '  Alpha-   =  ', a_aminus
222         WRITE(inum,*)
223      ENDIF
224
225 9006 FORMAT(1X,A,ES24.16)
226 9010 FORMAT(1X,A,ES12.5,A,F10.5,A)
227 9020 FORMAT(1X,A,F10.5,A)
228 9030 FORMAT(1X,A,F8.2,A)
229
230   END SUBROUTINE flx_fwb
231
232
233   SUBROUTINE flx_fwb_init
234      !!---------------------------------------------------------------------
235      !!                ***  ROUTINE flx_fwb_init  ***
236      !!
237      !! ** Purpose :
238      !!
239      !! ** Method  :
240      !! 
241      !! History :
242      !!   9.0  !  03-09  (G. Madec)  Original code
243      !!----------------------------------------------------------------------
244      !! * Local declarations
245      LOGICAL ::   llbon
246      CHARACTER (len=32) ::   &
247         clname = 'EMPave_old.dat'
248      INTEGER ::   inum = 11         ! temporary logical unit
249      INTEGER ::   iyear
250
251      NAMELIST/namfwb/ ln_fwb 
252      !!----------------------------------------------------------------------
253         
254      ! Read Namelist namfwb : freshWater Budget correction
255      ! --------------------
256      REWIND( numnam )
257      READ  ( numnam, namfwb )           
258     
259      ! Parameter control and print
260      ! ---------------------------
261      ! Control print
262      IF(lwp) THEN
263         WRITE(numout,*)
264         WRITE(numout,*) 'flx_fwb_init : FreshWater Budget correction'
265         WRITE(numout,*) '~~~~~~~~~~~~'
266         WRITE(numout,*) '               Namelist namfwb : set fwb parameters'
267         WRITE(numout,*) '                  use or not fwb correction      ln_fwb   = ', ln_fwb
268      ENDIF
269      ! Option consistency
270#if ! defined key_dynspg_fsc
271      IF(lwp) WRITE '               Rigid-lid option, fwb correction is useless, but valid'
272#endif
273      IF( lk_cpl ) THEN
274         IF(lwp) WRITE(numout,*) '               Coupled option, fwb correction is a flux correction ! '
275         IF(lwp) WRITE(numout,*) '               ln_fwb = .FALSE. is recommanded'
276      ENDIF
277
278      !                                        ! ==============================
279      IF( ln_fwb ) THEN                        !  Freshwater budget correction
280         !                                     ! ==============================
281         ! Read the corrective factor on precipitations (empold)
282         INQUIRE( FILE=clname, EXIST=llbon )
283         IF( llbon ) THEN
284            OPEN ( inum, FILE=clname)
285            READ ( inum, "(24X,I8,2ES24.16)" ) iyear, a_fwb_b, a_fwb
286            CLOSE( inum )
287            empold = a_fwb                  ! current year freshwater budget correction
288            !                               ! estimate from the previous year budget
289            IF(lwp)WRITE(numout,*)
290            IF(lwp)WRITE(numout,*)'flx_fwb_init : year = ',iyear  , ' freshwater budget correction = ', empold
291            IF(lwp)WRITE(numout,*)'               year = ',iyear-1, ' freshwater budget read       = ', a_fwb
292            IF(lwp)WRITE(numout,*)'               year = ',iyear-2, ' freshwater budget read       = ', a_fwb_b
293         ELSE
294            IF(lwp) WRITE(numout,*)
295            IF(lwp) WRITE(numout,*)'flx_fwb_init : unable to read the file', clname
296            nstop = nstop + 1
297         ENDIF
298         !                                    ! ==============================
299      ELSE                                    !      NO  budget correction
300         !                                    ! ==============================
301         IF(lwp) WRITE(numout,*)
302         IF(lwp) WRITE(numout,*)'flx_fwb_init : NO freshwater budget correction'
303         empold = 0.e0
304      ENDIF
305
306   END SUBROUTINE flx_fwb_init
307
308#else
309   !!----------------------------------------------------------------------
310   !!   Default case :                       
311   !!----------------------------------------------------------------------
312   LOGICAL, PUBLIC ::   ln_fwb = .FALSE.   !: no fwb forced
313CONTAINS
314   SUBROUTINE flx_fwb( kt )                ! dummy routine
315      WRITE(*,*) 'flx_fwb: You should not have seen this print! error?', kt
316   END SUBROUTINE flx_fwb
317   SUBROUTINE flx_fwb_init
318   END SUBROUTINE flx_fwb_init
319#endif
320   !!======================================================================
321END MODULE flxfwb
Note: See TracBrowser for help on using the repository browser.