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 branches/dev_002_LIM/NEMO/OPA_SRC/SBC – NEMO

source: branches/dev_002_LIM/NEMO/OPA_SRC/SBC/flxfwb.F90 @ 824

Last change on this file since 824 was 824, checked in by rblod, 16 years ago

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