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

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

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