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.
diafwb.F90 in branches/DEV_R1821_Rivers/NEMO/OPA_SRC/DIA – NEMO

source: branches/DEV_R1821_Rivers/NEMO/OPA_SRC/DIA/diafwb.F90 @ 1938

Last change on this file since 1938 was 1938, checked in by rfurner, 14 years ago

rnf has been separated from emp and emps. Also temperature and salinity of runoff can be specified, and runoff can be added to a user specified depth
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 18.4 KB
Line 
1MODULE diafwb
2   !!======================================================================
3   !!                       ***  MODULE  diafwb  ***
4   !! Ocean diagnostics: freshwater budget
5   !!======================================================================
6   !! History :  8.2  !  01-02  (E. Durand)  Original code
7   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
8   !!            9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
9   !!----------------------------------------------------------------------
10#if ( defined key_orca_r2 || defined  key_orca_r4 )  && ! defined key_coupled
11   !!----------------------------------------------------------------------
12   !!   "key_orca_r2 or 4"
13   !!----------------------------------------------------------------------
14   !!----------------------------------------------------------------------
15   !!   dia_fwb     : freshwater budget for global ocean configurations
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers
18   USE dom_oce         ! ocean space and time domain
19   USE phycst          ! physical constants
20   USE sbc_oce         ! ???
21   USE zdf_oce         ! ocean vertical physics
22   USE in_out_manager  ! I/O manager
23   USE lib_mpp         ! distributed memory computing library
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC dia_fwb    ! routine called by step.F90
29
30   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .TRUE.    !: fresh water budget flag
31
32   REAL(wp)               ::   a_fwf ,          &
33      &                        a_sshb, a_sshn, a_salb, a_saln
34   REAL(wp), DIMENSION(4) ::   a_flxi, a_flxo, a_temi, a_temo, a_sali, a_salo
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !!   OPA 9.0 , LOCEAN-IPSL (2006)
41   !! $Id$
42   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44
45CONTAINS
46
47   SUBROUTINE dia_fwb( kt )
48      !!---------------------------------------------------------------------
49      !!                  ***  ROUTINE dia_fwb  ***
50      !!     
51      !! ** Purpose :
52      !!----------------------------------------------------------------------
53      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
54      !!
55      INTEGER :: inum             ! temporary logical unit
56      INTEGER :: ji, jj, jk, jt   ! dummy loop indices
57      INTEGER :: ii0, ii1, ij0, ij1
58      REAL(wp) ::   zarea, zvol, zwei
59      REAL(wp) ::  ztemi(4), ztemo(4), zsali(4), zsalo(4), zflxi(4), zflxo(4)
60      REAL(wp) ::  zt, zs, zu 
61      REAL(wp) ::  zsm0, zfwfnew
62      !!----------------------------------------------------------------------
63
64      ! Mean global salinity
65      zsm0 = 34.72654
66
67      ! To compute fwf mean value mean fwf
68
69      IF( kt == nit000 ) THEN
70
71         a_fwf    = 0.e0
72         a_sshb   = 0.e0 ! valeur de ssh au debut de la simulation
73         a_salb   = 0.e0 ! valeur de sal au debut de la simulation
74         ! sshb used because diafwb called after tranxt (i.e. after the swap)
75         a_sshb = SUM( e1t(:,:) * e2t(:,:) * sshb(:,:) * tmask_i(:,:) )
76         IF( lk_mpp )   CALL mpp_sum( a_sshb )      ! sum over the global domain
77
78         DO jk = 1, jpkm1
79            DO jj = 2, jpjm1
80               DO ji = fs_2, fs_jpim1   ! vector opt.
81                  zwei  = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
82                  a_salb = a_salb + ( sb(ji,jj,jk) - zsm0 ) * zwei
83               END DO
84            END DO
85         END DO
86         IF( lk_mpp )   CALL mpp_sum( a_salb )      ! sum over the global domain
87      ENDIF
88     
89      a_fwf    = SUM( e1t(:,:) * e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) ) 
90      IF( lk_mpp )   CALL mpp_sum( a_fwf    )       ! sum over the global domain
91
92      IF( kt == nitend ) THEN
93         a_sshn = 0.e0
94         a_saln = 0.e0
95         zarea = 0.e0
96         zvol  = 0.e0
97         zfwfnew = 0.e0
98         ! Mean sea level at nitend
99         a_sshn = SUM( e1t(:,:) * e2t(:,:) * sshn(:,:) * tmask_i(:,:) )
100         IF( lk_mpp )   CALL mpp_sum( a_sshn )      ! sum over the global domain
101         zarea  = SUM( e1t(:,:) * e2t(:,:) *             tmask_i(:,:) )
102         IF( lk_mpp )   CALL mpp_sum( zarea  )      ! sum over the global domain
103         
104         DO jk = 1, jpkm1   
105            DO jj = 2, jpjm1
106               DO ji = fs_2, fs_jpim1   ! vector opt.
107                  zwei  = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
108                  a_saln = a_saln + ( sn(ji,jj,jk) - zsm0 ) * zwei
109                  zvol  = zvol  + zwei
110               END DO
111            END DO
112         END DO
113         IF( lk_mpp )   CALL mpp_sum( a_saln )      ! sum over the global domain
114         IF( lk_mpp )   CALL mpp_sum( zvol )      ! sum over the global domain
115         
116         ! Conversion in m3
117         a_fwf    = a_fwf * rdttra(1) * 1.e-3 
118         
119         ! fwf correction to bring back the mean ssh to zero
120         zfwfnew = a_sshn / ( ( nitend - nit000 + 1 ) * rdt ) * 1.e3 / zarea
121
122      ENDIF
123
124
125      ! Calcul des termes de transport
126      ! ------------------------------
127     
128      ! 1 --> Gibraltar
129      ! 2 --> Cadiz
130      ! 3 --> Red Sea
131      ! 4 --> Baltic Sea
132
133      IF( kt == nit000 ) THEN
134         a_flxi(:) = 0.e0
135         a_flxo(:) = 0.e0
136         a_temi(:) = 0.e0
137         a_temo(:) = 0.e0
138         a_sali(:) = 0.e0
139         a_salo(:) = 0.e0
140      ENDIF
141
142      zflxi(:) = 0.e0
143      zflxo(:) = 0.e0
144      ztemi(:) = 0.e0
145      ztemo(:) = 0.e0
146      zsali(:) = 0.e0
147      zsalo(:) = 0.e0
148
149      ! Mean flow at Gibraltar
150
151      IF( cp_cfg == "orca" ) THEN 
152               
153         SELECT CASE ( jp_cfg )
154         !                                           ! =======================
155         CASE ( 4 )                                  !  ORCA_R4 configuration
156            !                                        ! =======================
157            ii0 = 70   ;   ii1 = 70
158            ij0 = 52   ;   ij1 = 52
159            !                                        ! =======================
160         CASE ( 2 )                                  !  ORCA_R2 configuration
161            !                                        ! =======================
162            ii0 = 140   ;   ii1 = 140
163            ij0 = 102   ;   ij1 = 102
164            !                                        ! =======================
165         CASE DEFAULT                                !    ORCA R05 or R025
166            !                                        ! =======================
167            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
168            !
169         END SELECT
170         !
171         DO ji = mi0(ii0), mi1(ii1)
172            DO jj = mj0(ij0), mj1(ij1)
173               DO jk = 1, jpk 
174                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
175                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
176                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
177
178                  IF( un(ji,jj,jk) > 0.e0 ) THEN
179                     zflxi(1) = zflxi(1) +    zu
180                     ztemi(1) = ztemi(1) + zt*zu
181                     zsali(1) = zsali(1) + zs*zu
182                  ELSE
183                     zflxo(1) = zflxo(1) +    zu
184                     ztemo(1) = ztemo(1) + zt*zu
185                     zsalo(1) = zsalo(1) + zs*zu
186                  ENDIF
187               END DO
188            END DO
189         END DO
190      ENDIF
191     
192      ! Mean flow at Cadiz
193      IF( cp_cfg == "orca" ) THEN
194               
195         SELECT CASE ( jp_cfg )
196         !                                           ! =======================
197         CASE ( 4 )                                  !  ORCA_R4 configuration
198            !                                        ! =======================
199            ii0 = 69   ;   ii1 = 69
200            ij0 = 52   ;   ij1 = 52
201            !                                        ! =======================
202         CASE ( 2 )                                  !  ORCA_R2 configuration
203            !                                        ! =======================
204            ii0 = 137   ;   ii1 = 137
205            ij0 = 101   ;   ij1 = 102
206            !                                        ! =======================
207         CASE DEFAULT                                !    ORCA R05 or R025
208            !                                        ! =======================
209            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
210            !
211         END SELECT
212         !
213         DO ji = mi0(ii0), mi1(ii1)
214            DO jj = mj0(ij0), mj1(ij1)
215               DO jk = 1, jpk 
216                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
217                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
218                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
219                 
220                  IF( un(ji,jj,jk) > 0.e0 ) THEN
221                     zflxi(2) = zflxi(2) +    zu
222                     ztemi(2) = ztemi(2) + zt*zu
223                     zsali(2) = zsali(2) + zs*zu
224                  ELSE
225                     zflxo(2) = zflxo(2) +    zu
226                     ztemo(2) = ztemo(2) + zt*zu
227                     zsalo(2) = zsalo(2) + zs*zu
228                  ENDIF
229               END DO
230            END DO
231         END DO
232      ENDIF
233
234      ! Mean flow at Red Sea entrance
235      IF( cp_cfg == "orca" ) THEN
236               
237         SELECT CASE ( jp_cfg )
238         !                                           ! =======================
239         CASE ( 4 )                                  !  ORCA_R4 configuration
240            !                                        ! =======================
241            ii0 = 83   ;   ii1 = 83
242            ij0 = 45   ;   ij1 = 45
243            !                                        ! =======================
244         CASE ( 2 )                                  !  ORCA_R2 configuration
245            !                                        ! =======================
246            ii0 = 160   ;   ii1 = 160
247            ij0 = 88    ;   ij1 = 88 
248            !                                        ! =======================
249         CASE DEFAULT                                !    ORCA R05 or R025
250            !                                        ! =======================
251            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
252            !
253         END SELECT
254         !
255         DO ji = mi0(ii0), mi1(ii1)
256            DO jj = mj0(ij0), mj1(ij1)
257               DO jk = 1, jpk 
258                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
259                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
260                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
261                 
262                  IF( un(ji,jj,jk) > 0.e0 ) THEN
263                     zflxi(3) = zflxi(3) +    zu
264                     ztemi(3) = ztemi(3) + zt*zu
265                     zsali(3) = zsali(3) + zs*zu
266                  ELSE
267                     zflxo(3) = zflxo(3) +    zu
268                     ztemo(3) = ztemo(3) + zt*zu
269                     zsalo(3) = zsalo(3) + zs*zu
270                  ENDIF
271               END DO
272            END DO
273         END DO
274      ENDIF
275
276      ! Mean flow at Baltic Sea entrance
277      IF( cp_cfg == "orca" ) THEN
278               
279         SELECT CASE ( jp_cfg )
280         !                                           ! =======================
281         CASE ( 4 )                                  !  ORCA_R4 configuration
282            !                                        ! =======================
283            ii0 = 1     ;   ii1 = 1 
284            ij0 = 1     ;   ij1 = 1 
285            !                                        ! =======================
286         CASE ( 2 )                                  !  ORCA_R2 configuration
287            !                                        ! =======================
288            ii0 = 146   ;   ii1 = 146 
289            ij0 = 116   ;   ij1 = 116
290            !                                        ! =======================
291         CASE DEFAULT                                !    ORCA R05 or R025
292            !                                        ! =======================
293            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
294            !
295         END SELECT
296         !
297         DO ji = mi0(ii0), mi1(ii1)
298            DO jj = mj0(ij0), mj1(ij1)
299               DO jk = 1, jpk
300                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
301                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
302                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
303                 
304                  IF( un(ji,jj,jk) > 0.e0 ) THEN
305                     zflxi(4) = zflxi(4) +    zu
306                     ztemi(4) = ztemi(4) + zt*zu
307                     zsali(4) = zsali(4) + zs*zu
308                  ELSE
309                     zflxo(4) = zflxo(4) +    zu
310                     ztemo(4) = ztemo(4) + zt*zu
311                     zsalo(4) = zsalo(4) + zs*zu
312                  ENDIF
313               END DO
314            END DO
315         END DO
316      ENDIF
317
318      ! Sum at each time-step
319      DO jt = 1, 4 
320         !
321         IF( zflxi(jt) /= 0.e0 ) THEN
322            a_flxi(jt) = a_flxi(jt) + zflxi(jt)
323            a_temi(jt) = a_temi(jt) + ztemi(jt)/zflxi(jt)
324            a_sali(jt) = a_sali(jt) + zsali(jt)/zflxi(jt)
325         ENDIF
326         !
327         IF( zflxo(jt) /= 0.e0 ) THEN
328            a_flxo(jt) = a_flxo(jt) + zflxo(jt)
329            a_temo(jt) = a_temo(jt) + ztemo(jt)/zflxo(jt)
330            a_salo(jt) = a_salo(jt) + zsalo(jt)/zflxo(jt)
331         ENDIF
332         !
333      END DO
334
335      IF( kt == nitend ) THEN
336         DO jt = 1, 4 
337            a_flxi(jt) = a_flxi(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
338            a_temi(jt) = a_temi(jt) /   FLOAT( nitend - nit000 + 1 )
339            a_sali(jt) = a_sali(jt) /   FLOAT( nitend - nit000 + 1 )
340            a_flxo(jt) = a_flxo(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
341            a_temo(jt) = a_temo(jt) /   FLOAT( nitend - nit000 + 1 )
342            a_salo(jt) = a_salo(jt) /   FLOAT( nitend - nit000 + 1 )
343         END DO
344         IF( lk_mpp ) THEN
345            CALL mpp_sum( a_flxi, 4 )      ! sum over the global domain
346            CALL mpp_sum( a_temi, 4 )      ! sum over the global domain
347            CALL mpp_sum( a_sali, 4 )      ! sum over the global domain
348
349            CALL mpp_sum( a_flxo, 4 )      ! sum over the global domain
350            CALL mpp_sum( a_temo, 4 )      ! sum over the global domain
351            CALL mpp_sum( a_salo, 4 )      ! sum over the global domain
352         ENDIF
353      ENDIF
354
355
356      ! Ecriture des diagnostiques
357      ! --------------------------
358
359      IF ( kt == nitend .AND. cp_cfg == "orca" .AND. lwp ) THEN
360
361         CALL ctl_opn( inum, 'STRAIT.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
362         WRITE(inum,*)
363         WRITE(inum,*)    'Net freshwater budget '
364         WRITE(inum,9010) '  fwf    = ',a_fwf,   ' m3 =', a_fwf   /(FLOAT(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
365         WRITE(inum,*)
366         WRITE(inum,9010) '  zarea =',zarea
367         WRITE(inum,9010) '  zvol  =',zvol
368         WRITE(inum,*)
369         WRITE(inum,*)    'Mean sea level : '
370         WRITE(inum,9010) '  at nit000 = ',a_sshb        ,' m3 '
371         WRITE(inum,9010) '  at nitend = ',a_sshn        ,' m3 '
372         WRITE(inum,9010) '  diff      = ',(a_sshn-a_sshb),' m3 =', (a_sshn-a_sshb)/(FLOAT(nitend-nit000+1)*rdt) * 1.e-6,' Sv'
373         WRITE(inum,9020) '  mean sea level elevation    =', a_sshn/zarea,' m'
374         WRITE(inum,*)
375         WRITE(inum,*)    'Anomaly of salinity content : '
376         WRITE(inum,9010) '  at nit000 = ',a_salb        ,' psu.m3 '
377         WRITE(inum,9010) '  at nitend = ',a_saln        ,' psu.m3 '
378         WRITE(inum,9010) '  diff      = ',(a_saln-a_salb),' psu.m3'
379         WRITE(inum,*)
380         WRITE(inum,*)    'Mean salinity : '
381         WRITE(inum,9020) '  at nit000 =',a_salb/zvol+zsm0   ,' psu '
382         WRITE(inum,9020) '  at nitend =',a_saln/zvol+zsm0   ,' psu '
383         WRITE(inum,9020) '  diff      =',(a_saln-a_salb)/zvol,' psu'
384         WRITE(inum,9020) '  S-SLevitus=',a_saln/zvol,' psu'
385         WRITE(inum,*)
386         WRITE(inum,*)    'Gibraltar : '
387         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(1)
388         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(1)
389         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(1)
390         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(1)
391         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(1)
392         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(1)
393         WRITE(inum,*)
394         WRITE(inum,*)    'Cadiz : '
395         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(2)
396         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(2)
397         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(2)
398         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(2)
399         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(2)
400         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(2)
401         WRITE(inum,*)
402         WRITE(inum,*)    'Bab el Mandeb : '
403         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(3)
404         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(3)
405         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(3)
406         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(3)
407         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(3)
408         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(3)
409         WRITE(inum,*)
410         WRITE(inum,*)    'Baltic : '
411         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(4)
412         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(4)
413         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(4)
414         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(4)
415         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(4)
416         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(4)
417         CLOSE(inum)
418      ENDIF
419
420 9005 FORMAT(1X,A,ES24.16)
421 9010 FORMAT(1X,A,ES12.5,A,F10.5,A)
422 9020 FORMAT(1X,A,F10.5,A)
423 9030 FORMAT(1X,A,F9.4,A)
424
425   END SUBROUTINE dia_fwb
426
427#else
428   !!----------------------------------------------------------------------
429   !!   Default option :                                       Dummy Module
430   !!----------------------------------------------------------------------
431   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .FALSE.    !: fresh water budget flag
432CONTAINS
433   SUBROUTINE dia_fwb( kt )        ! Empty routine
434      WRITE(*,*) 'dia_fwb: : You should not have seen this print! error?', kt
435   END SUBROUTINE dia_fwb
436#endif
437
438   !!======================================================================
439END MODULE diafwb
Note: See TracBrowser for help on using the repository browser.