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 on Ticket #127 – Attachment – NEMO

Ticket #127: diafwb.F90

File diafwb.F90, 27.0 KB (added by nemo_user, 16 years ago)

Modified diafwb.F90 source file

Line 
1MODULE diafwb
2   !!======================================================================
3   !!                       ***  MODULE  diafwb  ***
4   !! Ocean diagnostics: freshwater budget
5   !!======================================================================
6#if ( defined key_orca_r2 || defined  key_orca_r4 ) && ! defined key_dynspg_rl && ! defined key_coupled
7   !!----------------------------------------------------------------------
8   !!   NOT "key_dynspg_rl" and "key_orca_r2 or 4"
9   !!----------------------------------------------------------------------
10   !!   dia_fwb     : 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 zdf_oce         ! ocean vertical physics
17   USE in_out_manager  ! I/O manager
18   USE flxrnf          ! ???
19   USE ocesbc          ! ???
20   USE blk_oce         ! ???
21   USE flxblk          ! atmospheric surface quantity
22   USE lib_mpp         ! distributed memory computing library
23
24   IMPLICIT NONE
25   PRIVATE
26
27   !! * Routine accessibility
28   PUBLIC dia_fwb    ! routine called by step.F90
29
30   !! * Shared module variables
31   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .TRUE.    !: fresh water budget flag
32
33   !! * Module variables
34   REAL(wp) ::   &
35      a_emp , a_precip, a_rnf,   &
36      a_sshb, a_sshn, a_salb, a_saln,   &
37      a_aminus, a_aplus
38   REAL(wp), DIMENSION(4) ::   &
39      a_flxi, a_flxo, a_temi, a_temo, a_sali, a_salo
40
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !!   OPA 9.0 , LOCEAN-IPSL (2005)
46   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/DIA/diafwb.F90,v 1.11 2007/06/29 17:01:51 opalod Exp $
47   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
48   !!----------------------------------------------------------------------
49
50CONTAINS
51
52   SUBROUTINE dia_fwb( kt )
53      !!---------------------------------------------------------------------
54      !!                  ***  ROUTINE dia_fwb  ***
55      !!     
56      !! ** Purpose :
57      !!
58      !! ** Method :
59      !!
60      !! History :
61      !!   8.2  !  01-02  (E. Durand)  Original code
62      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
63      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
64      !!----------------------------------------------------------------------
65      !! * Arguments
66      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
67
68      !! * Local declarations
69      CHARACTER (len=32) ::  clname
70      INTEGER :: inum             ! temporary logical unit
71      INTEGER :: ji, jj, jk, jt   ! dummy loop indices
72      INTEGER :: ii0, ii1, ij0, ij1
73      INTEGER :: ilimit ! Adjustable loop limit to allow for edge cases
74      REAL(wp) ::   zarea, zvol, zwei
75      REAL(wp) ::  ztemi(4), ztemo(4), zsali(4), zsalo(4), zflxi(4), zflxo(4)
76      REAL(wp) ::  zt, zs, zu 
77      REAL(wp) ::  zsm0, zempnew
78      !!----------------------------------------------------------------------
79
80      ! Mean global salinity
81      zsm0 = 34.72654
82
83      ! To compute emp mean value mean emp
84
85      IF( kt == nit000 ) THEN
86
87         a_emp    = 0.e0
88         a_precip = 0.e0
89         a_rnf    = 0.e0
90         a_sshb   = 0.e0 ! valeur de ssh au debut de la simulation
91         a_salb   = 0.e0 ! valeur de sal au debut de la simulation
92         a_aminus = 0.e0
93         a_aplus  = 0.e0
94         ! sshb used because diafwb called after tranxt (i.e. after the swap)
95         a_sshb = SUM( e1t(:,:) * e2t(:,:) * sshb(:,:) * tmask_i(:,:) )
96         IF( lk_mpp )   CALL mpp_sum( a_sshb )      ! sum over the global domain
97
98         DO jk = 1, jpkm1
99            DO jj = 2, jpjm1
100               DO ji = fs_2, fs_jpim1   ! vector opt.
101                  zwei  = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
102                  a_salb = a_salb + ( sb(ji,jj,jk) - zsm0 ) * zwei
103               END DO
104            END DO
105         END DO
106         IF( lk_mpp )   CALL mpp_sum( a_salb )      ! sum over the global domain
107      ENDIF
108     
109      a_emp    = SUM( e1t(:,:) * e2t(:,:) * emp   (:,:) * tmask_i(:,:) )
110      IF( lk_mpp )   CALL mpp_sum( a_emp    )       ! sum over the global domain
111#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
112      a_precip = SUM( e1t(:,:) * e2t(:,:) * watm  (:,:) * tmask_i(:,:) )
113      IF( lk_mpp )   CALL mpp_sum( a_precip )       ! sum over the global domain
114#endif
115      a_rnf    = SUM( e1t(:,:) * e2t(:,:) * runoff(:,:) * tmask_i(:,:) )
116      IF( lk_mpp )   CALL mpp_sum( a_rnf    )       ! sum over the global domain
117
118      IF( aminus /= 0.e0 ) a_aminus = a_aminus + ( MIN( aplus, aminus ) / aminus )
119      IF( aplus  /= 0.e0 ) a_aplus  = a_aplus  + ( MIN( aplus, aminus ) / aplus  )
120
121      IF( kt == nitend ) THEN
122         a_sshn = 0.e0
123         a_saln = 0.e0
124         zarea = 0.e0
125         zvol  = 0.e0
126         zempnew = 0.e0
127         ! Mean sea level at nitend
128         a_sshn = SUM( e1t(:,:) * e2t(:,:) * sshn(:,:) * tmask_i(:,:) )
129         IF( lk_mpp )   CALL mpp_sum( a_sshn )      ! sum over the global domain
130         zarea  = SUM( e1t(:,:) * e2t(:,:) *             tmask_i(:,:) )
131         IF( lk_mpp )   CALL mpp_sum( zarea  )      ! sum over the global domain
132         
133         DO jk = 1, jpkm1   
134            DO jj = 2, jpjm1
135               DO ji = fs_2, fs_jpim1   ! vector opt.
136                  zwei  = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
137                  a_saln = a_saln + ( sn(ji,jj,jk) - zsm0 ) * zwei
138                  zvol  = zvol  + zwei
139               END DO
140            END DO
141         END DO
142         IF( lk_mpp )   CALL mpp_sum( a_saln )      ! sum over the global domain
143         IF( lk_mpp )   CALL mpp_sum( zvol )      ! sum over the global domain
144         
145         a_aminus = a_aminus / ( nitend - nit000 + 1 )
146         a_aplus  = a_aplus  / ( nitend - nit000 + 1 )
147
148         ! Conversion in m3
149         a_emp    = a_emp * rdttra(1) * 1.e-3 
150         a_precip = a_precip * rdttra(1) * 1.e-3 / rday
151         a_rnf    = a_rnf * rdttra(1) * 1.e-3
152         
153         ! Alpha1=Alpha0-Rest/(Precip+runoff)
154         !  C A U T I O N : precipitations are negative !!     
155         
156         zempnew = a_sshn / ( ( nitend - nit000 + 1 ) * rdt ) * 1.e3 / zarea
157
158      ENDIF
159
160
161      ! Calcul des termes de transport
162      ! ------------------------------
163     
164      ! 1 --> Gibraltar
165      ! 2 --> Cadiz
166      ! 3 --> Red Sea
167      ! 4 --> Baltic Sea
168
169      IF( kt == nit000 ) THEN
170         a_flxi(:) = 0.e0
171         a_flxo(:) = 0.e0
172         a_temi(:) = 0.e0
173         a_temo(:) = 0.e0
174         a_sali(:) = 0.e0
175         a_salo(:) = 0.e0
176      ENDIF
177
178      zflxi(:) = 0.e0
179      zflxo(:) = 0.e0
180      ztemi(:) = 0.e0
181      ztemo(:) = 0.e0
182      zsali(:) = 0.e0
183      zsalo(:) = 0.e0
184
185      ! Mean flow at Gibraltar
186
187      IF( cp_cfg == "orca" ) THEN 
188               
189         SELECT CASE ( jp_cfg )
190         !                                           ! =======================
191         CASE ( 4 )                                  !  ORCA_R4 configuration
192            !                                        ! =======================
193            ii0 = 70   ;   ii1 = 70
194            ij0 = 52   ;   ij1 = 52
195            !                                        ! =======================
196         CASE ( 2 )                                  !  ORCA_R2 configuration
197            !                                        ! =======================
198            ii0 = 139   ;   ii1 = 139
199            ij0 = 102   ;   ij1 = 102
200            !                                        ! =======================
201         CASE DEFAULT                                !    ORCA R05 or R025
202            !                                        ! =======================
203            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
204            !
205         END SELECT
206         !
207!!$         DO ji = mi0(ii0), mi1(ii1)
208!!$            DO jj = mj0(ij0), mj1(ij1)
209!!$               DO jk = 1, 18
210!!$                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
211!!$                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
212!!$                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
213!!$
214!!$                  IF( un(ji,jj,jk) > 0.e0 ) THEN
215!!$                     zflxi(1) = zflxi(1) +    zu
216!!$                     ztemi(1) = ztemi(1) + zt*zu
217!!$                     zsali(1) = zsali(1) + zs*zu
218!!$                  ELSE
219!!$                     zflxo(1) = zflxo(1) +    zu
220!!$                     ztemo(1) = ztemo(1) + zt*zu
221!!$                     zsalo(1) = zsalo(1) + zs*zu
222!!$                  ENDIF
223!!$               END DO
224!!$            END DO
225!!$         END DO
226
227         ! ARPDBG Correct loop limits if at edge
228         IF(mi1(ii1) == jpi)THEN
229            ilimit = jpi - 1
230         ELSE
231            ilimit = mi1(ii1)
232         END IF
233
234         DO ji = mi0(ii0), ilimit
235            DO jj = mj0(ij0), mj1(ij1)
236               DO jk = 1, 18 
237                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
238                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
239                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
240
241                  IF( un(ji,jj,jk) > 0.e0 ) THEN
242                     zflxi(1) = zflxi(1) +    zu
243                     ztemi(1) = ztemi(1) + zt*zu
244                     zsali(1) = zsali(1) + zs*zu
245                  ELSE
246                     zflxo(1) = zflxo(1) +    zu
247                     ztemo(1) = ztemo(1) + zt*zu
248                     zsalo(1) = zsalo(1) + zs*zu
249                  ENDIF
250               END DO
251            END DO
252         END DO
253
254         ! ARPDBG - deal with edge case separately
255         IF(mi1(ii1) == jpi)THEN
256
257            ji = jpi
258            DO jk = 1, 18 
259               DO jj = mj0(ij0), mj1(ij1)
260                  zt = tn(ji,jj,jk)
261                  zs = sn(ji,jj,jk)
262                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
263
264                  IF( un(ji,jj,jk) > 0.e0 ) THEN
265                     zflxi(1) = zflxi(1) +    zu
266                     ztemi(1) = ztemi(1) + zt*zu
267                     zsali(1) = zsali(1) + zs*zu
268                  ELSE
269                     zflxo(1) = zflxo(1) +    zu
270                     ztemo(1) = ztemo(1) + zt*zu
271                     zsalo(1) = zsalo(1) + zs*zu
272                  ENDIF
273               END DO
274            END DO
275
276         ENDIF ! Edge case
277
278      ENDIF
279     
280      ! Mean flow at Cadiz
281      IF( cp_cfg == "orca" ) THEN
282               
283         SELECT CASE ( jp_cfg )
284         !                                           ! =======================
285         CASE ( 4 )                                  !  ORCA_R4 configuration
286            !                                        ! =======================
287            ii0 = 69   ;   ii1 = 69
288            ij0 = 52   ;   ij1 = 52
289            !                                        ! =======================
290         CASE ( 2 )                                  !  ORCA_R2 configuration
291            !                                        ! =======================
292            ii0 = 137   ;   ii1 = 137
293            ij0 = 102   ;   ij1 = 102
294            !                                        ! =======================
295         CASE DEFAULT                                !    ORCA R05 or R025
296            !                                        ! =======================
297            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
298            !
299         END SELECT
300         !
301!!$         DO ji = mi0(ii0), mi1(ii1)
302!!$            DO jj = mj0(ij0), mj1(ij1)
303!!$               DO jk = 1, 23
304!!$                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
305!!$                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
306!!$                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
307!!$                 
308!!$                  IF( un(ji,jj,jk) > 0.e0 ) THEN
309!!$                     zflxi(2) = zflxi(2) +    zu
310!!$                     ztemi(2) = ztemi(2) + zt*zu
311!!$                     zsali(2) = zsali(2) + zs*zu
312!!$                  ELSE
313!!$                     zflxo(2) = zflxo(2) +    zu
314!!$                     ztemo(2) = ztemo(2) + zt*zu
315!!$                     zsalo(2) = zsalo(2) + zs*zu
316!!$                  ENDIF
317!!$               END DO
318!!$            END DO
319!!$         END DO
320
321         ! ARPDBG Correct loop limits if at edge
322         IF(mi1(ii1) == jpi)THEN
323            ilimit = jpi - 1
324         ELSE
325            ilimit = mi1(ii1)
326         END IF
327
328         DO jk = 1, 23 
329            DO jj = mj0(ij0), mj1(ij1)
330               DO ji = mi0(ii0), ilimit
331                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
332                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
333                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
334                 
335                  IF( un(ji,jj,jk) > 0.e0 ) THEN
336                     zflxi(2) = zflxi(2) +    zu
337                     ztemi(2) = ztemi(2) + zt*zu
338                     zsali(2) = zsali(2) + zs*zu
339                  ELSE
340                     zflxo(2) = zflxo(2) +    zu
341                     ztemo(2) = ztemo(2) + zt*zu
342                     zsalo(2) = zsalo(2) + zs*zu
343                  ENDIF
344               END DO
345            END DO
346         END DO
347
348         ! ARPDBG - deal with edge case separately
349         IF(mi1(ii1) == jpi)THEN
350            ji = jpi
351
352            DO jk = 1, 23 
353               DO jj = mj0(ij0), mj1(ij1)
354
355                  zt = tn(ji,jj,jk)
356                  zs = sn(ji,jj,jk)
357                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
358                 
359                  IF( un(ji,jj,jk) > 0.e0 ) THEN
360                     zflxi(2) = zflxi(2) +    zu
361                     ztemi(2) = ztemi(2) + zt*zu
362                     zsali(2) = zsali(2) + zs*zu
363                  ELSE
364                     zflxo(2) = zflxo(2) +    zu
365                     ztemo(2) = ztemo(2) + zt*zu
366                     zsalo(2) = zsalo(2) + zs*zu
367                  ENDIF
368               END DO
369            END DO
370         END IF ! Edge case
371
372      ENDIF ! Cadiz
373
374      ! Mean flow at Red Sea entrance
375      IF( cp_cfg == "orca" ) THEN
376               
377         SELECT CASE ( jp_cfg )
378         !                                           ! =======================
379         CASE ( 4 )                                  !  ORCA_R4 configuration
380            !                                        ! =======================
381            ii0 = 83   ;   ii1 = 83
382            ij0 = 45   ;   ij1 = 45
383            !                                        ! =======================
384         CASE ( 2 )                                  !  ORCA_R2 configuration
385            !                                        ! =======================
386            ii0 = 161   ;   ii1 = 161
387            ij0 = 88    ;   ij1 = 88 
388            !                                        ! =======================
389         CASE DEFAULT                                !    ORCA R05 or R025
390            !                                        ! =======================
391            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
392            !
393         END SELECT
394         !
395!!$         DO ji = mi0(ii0), mi1(ii1)
396!!$            DO jj = mj0(ij0), mj1(ij1)
397!!$               DO jk = 1, 15
398!!$                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
399!!$                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
400!!$                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
401!!$                 
402!!$                  IF( un(ji,jj,jk) > 0.e0 ) THEN
403!!$                     zflxi(3) = zflxi(3) +    zu
404!!$                     ztemi(3) = ztemi(3) + zt*zu
405!!$                     zsali(3) = zsali(3) + zs*zu
406!!$                  ELSE
407!!$                     zflxo(3) = zflxo(3) +    zu
408!!$                     ztemo(3) = ztemo(3) + zt*zu
409!!$                     zsalo(3) = zsalo(3) + zs*zu
410!!$                  ENDIF
411!!$               END DO
412!!$            END DO
413!!$         END DO
414
415         ! ARPDBG Correct loop limits if at edge
416         IF(mi1(ii1) == jpi)THEN
417            ilimit = jpi - 1
418         ELSE
419            ilimit = mi1(ii1)
420         END IF
421
422         DO jk = 1, 15 
423            DO jj = mj0(ij0), mj1(ij1)
424               DO ji = mi0(ii0), ilimit
425                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
426                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
427                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
428                 
429                  IF( un(ji,jj,jk) > 0.e0 ) THEN
430                     zflxi(3) = zflxi(3) +    zu
431                     ztemi(3) = ztemi(3) + zt*zu
432                     zsali(3) = zsali(3) + zs*zu
433                  ELSE
434                     zflxo(3) = zflxo(3) +    zu
435                     ztemo(3) = ztemo(3) + zt*zu
436                     zsalo(3) = zsalo(3) + zs*zu
437                  ENDIF
438               END DO
439            END DO
440         END DO
441
442         ! ARPDBG - deal with edge case separately
443         IF(mi1(ii1) == jpi)THEN
444            ji = jpi
445            DO jk = 1, 15 
446               DO jj = mj0(ij0), mj1(ij1)
447                  zt = tn(ji,jj,jk)
448                  zs = sn(ji,jj,jk)
449                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
450                 
451                  IF( un(ji,jj,jk) > 0.e0 ) THEN
452                     zflxi(3) = zflxi(3) +    zu
453                     ztemi(3) = ztemi(3) + zt*zu
454                     zsali(3) = zsali(3) + zs*zu
455                  ELSE
456                     zflxo(3) = zflxo(3) +    zu
457                     ztemo(3) = ztemo(3) + zt*zu
458                     zsalo(3) = zsalo(3) + zs*zu
459                  ENDIF
460               END DO
461            END DO
462         END IF ! Edge case
463
464      ENDIF
465
466      ! Mean flow at Baltic Sea entrance
467      IF( cp_cfg == "orca" ) THEN
468               
469         SELECT CASE ( jp_cfg )
470         !                                           ! =======================
471         CASE ( 4 )                                  !  ORCA_R4 configuration
472            !                                        ! =======================
473            ii0 = 1     ;   ii1 = 1 
474            ij0 = 1     ;   ij1 = 1 
475            !                                        ! =======================
476         CASE ( 2 )                                  !  ORCA_R2 configuration
477            !                                        ! =======================
478            ii0 = 146   ;   ii1 = 146 
479            ij0 = 116   ;   ij1 = 116
480            !                                        ! =======================
481         CASE DEFAULT                                !    ORCA R05 or R025
482            !                                        ! =======================
483            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
484            !
485         END SELECT
486         !
487!!$         DO ji = mi0(ii0), mi1(ii1)
488!!$            DO jj = mj0(ij0), mj1(ij1)
489!!$               DO jk = 1, 20
490!!$                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
491!!$                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
492!!$                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
493!!$                 
494!!$                  IF( un(ji,jj,jk) > 0.e0 ) THEN
495!!$                     zflxi(4) = zflxi(4) +    zu
496!!$                     ztemi(4) = ztemi(4) + zt*zu
497!!$                     zsali(4) = zsali(4) + zs*zu
498!!$                  ELSE
499!!$                     zflxo(4) = zflxo(4) +    zu
500!!$                     ztemo(4) = ztemo(4) + zt*zu
501!!$                     zsalo(4) = zsalo(4) + zs*zu
502!!$                  ENDIF
503!!$               END DO
504!!$            END DO
505!!$         END DO
506
507         ! ARPDBG Correct loop limits if at edge
508         IF(mi1(ii1) == jpi)THEN
509            ilimit = jpi - 1
510         ELSE
511            ilimit = mi1(ii1)
512         END IF
513
514         DO jk = 1, 20
515            DO jj = mj0(ij0), mj1(ij1)
516               DO ji = mi0(ii0), ilimit
517                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
518                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
519                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
520                 
521                  IF( un(ji,jj,jk) > 0.e0 ) THEN
522                     zflxi(4) = zflxi(4) +    zu
523                     ztemi(4) = ztemi(4) + zt*zu
524                     zsali(4) = zsali(4) + zs*zu
525                  ELSE
526                     zflxo(4) = zflxo(4) +    zu
527                     ztemo(4) = ztemo(4) + zt*zu
528                     zsalo(4) = zsalo(4) + zs*zu
529                  ENDIF
530               END DO
531            END DO
532         END DO
533
534         ! ARPDBG - deal with edge case separately
535         IF(mi1(ii1) == jpi)THEN
536            ji = jpi
537
538            DO jk = 1, 20
539               DO jj = mj0(ij0), mj1(ij1)
540                  zt = tn(ji,jj,jk)
541                  zs = sn(ji,jj,jk)
542                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
543                 
544                  IF( un(ji,jj,jk) > 0.e0 ) THEN
545                     zflxi(4) = zflxi(4) +    zu
546                     ztemi(4) = ztemi(4) + zt*zu
547                     zsali(4) = zsali(4) + zs*zu
548                  ELSE
549                     zflxo(4) = zflxo(4) +    zu
550                     ztemo(4) = ztemo(4) + zt*zu
551                     zsalo(4) = zsalo(4) + zs*zu
552                  ENDIF
553               END DO
554            END DO
555         END IF ! Edge case
556
557      ENDIF
558
559      ! Sum at each time-step
560      DO jt = 1, 4 
561         IF( zflxi(jt) /= 0.e0 .AND. zflxo(jt) /= 0.e0 ) THEN
562            a_flxi(jt) = a_flxi(jt) + zflxi(jt)
563            a_temi(jt) = a_temi(jt) + ztemi(jt)/zflxi(jt)
564            a_sali(jt) = a_sali(jt) + zsali(jt)/zflxi(jt)
565            a_flxo(jt) = a_flxo(jt) + zflxo(jt)
566            a_temo(jt) = a_temo(jt) + ztemo(jt)/zflxo(jt)
567            a_salo(jt) = a_salo(jt) + zsalo(jt)/zflxo(jt)
568         ENDIF
569      END DO
570
571      IF( kt == nitend ) THEN
572         DO jt = 1, 4 
573            a_flxi(jt) = a_flxi(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
574            a_temi(jt) = a_temi(jt) /   FLOAT( nitend - nit000 + 1 )
575            a_sali(jt) = a_sali(jt) /   FLOAT( nitend - nit000 + 1 )
576            a_flxo(jt) = a_flxo(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
577            a_temo(jt) = a_temo(jt) /   FLOAT( nitend - nit000 + 1 )
578            a_salo(jt) = a_salo(jt) /   FLOAT( nitend - nit000 + 1 )
579         END DO
580         IF( lk_mpp ) THEN
581            CALL mpp_sum( a_flxi, 4 )      ! sum over the global domain
582            CALL mpp_sum( a_temi, 4 )      ! sum over the global domain
583            CALL mpp_sum( a_sali, 4 )      ! sum over the global domain
584
585            CALL mpp_sum( a_flxo, 4 )      ! sum over the global domain
586            CALL mpp_sum( a_temo, 4 )      ! sum over the global domain
587            CALL mpp_sum( a_salo, 4 )      ! sum over the global domain
588         ENDIF
589      ENDIF
590
591
592      ! Ecriture des diagnostiques
593      ! --------------------------
594
595      IF ( kt == nitend .AND. cp_cfg == "orca" ) THEN
596
597         clname = 'STRAIT.dat'
598         CALL ctlopn( inum, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
599            &         1, numout, lwp, 1 )
600         WRITE(inum,*)
601         WRITE(inum,*)    'Net freshwater budget '
602         WRITE(inum,9010) '  emp    = ',a_emp,   ' m3 =', a_emp   /(FLOAT(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
603         WRITE(inum,9010) '  precip = ',a_precip,' m3 =', a_precip/(FLOAT(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
604         WRITE(inum,9010) '  a_rnf  = ',a_rnf,   ' m3 =', a_rnf   /(FLOAT(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
605         WRITE(inum,*)
606         WRITE(inum,9010) '  zarea =',zarea
607         WRITE(inum,9010) '  zvol  =',zvol
608         WRITE(inum,*)
609         WRITE(inum,*)    'Mean sea level : '
610         WRITE(inum,9010) '  at nit000 = ',a_sshb        ,' m3 '
611         WRITE(inum,9010) '  at nitend = ',a_sshn        ,' m3 '
612         WRITE(inum,9010) '  diff      = ',(a_sshn-a_sshb),' m3 =', (a_sshn-a_sshb)/(FLOAT(nitend-nit000+1)*rdt) * 1.e-6,' Sv'
613         WRITE(inum,9020) '  mean sea level elevation    =', a_sshn/zarea,' m'
614         WRITE(inum,*)
615         WRITE(inum,*)    'Anomaly of salinity content : '
616         WRITE(inum,9010) '  at nit000 = ',a_salb        ,' psu.m3 '
617         WRITE(inum,9010) '  at nitend = ',a_saln        ,' psu.m3 '
618         WRITE(inum,9010) '  diff      = ',(a_saln-a_salb),' psu.m3'
619         WRITE(inum,*)
620         WRITE(inum,*)    'Mean salinity : '
621         WRITE(inum,9020) '  at nit000 =',a_salb/zvol+zsm0   ,' psu '
622         WRITE(inum,9020) '  at nitend =',a_saln/zvol+zsm0   ,' psu '
623         WRITE(inum,9020) '  diff      =',(a_saln-a_salb)/zvol,' psu'
624         WRITE(inum,9020) '  S-SLevitus=',a_saln/zvol,' psu'
625         WRITE(inum,*)
626         WRITE(inum,*)    'Coeff : '
627         WRITE(inum,9030) '  Alpha+   =  ', a_aplus
628         WRITE(inum,9030) '  Alpha-   =  ', a_aminus
629         WRITE(inum,*)
630         WRITE(inum,*)
631         WRITE(inum,*)    'Gibraltar : '
632         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(1)
633         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(1)
634         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(1)
635         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(1)
636         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(1)
637         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(1)
638         WRITE(inum,*)
639         WRITE(inum,*)    'Cadiz : '
640         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(2)
641         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(2)
642         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(2)
643         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(2)
644         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(2)
645         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(2)
646         WRITE(inum,*)
647         WRITE(inum,*)    'Bab el Mandeb : '
648         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(3)
649         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(3)
650         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(3)
651         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(3)
652         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(3)
653         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(3)
654         WRITE(inum,*)
655         WRITE(inum,*)    'Baltic : '
656         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(4)
657         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(4)
658         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(4)
659         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(4)
660         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(4)
661         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(4)
662         CLOSE(inum)
663      ENDIF
664
665 9005 FORMAT(1X,A,ES24.16)
666 9010 FORMAT(1X,A,ES12.5,A,F10.5,A)
667 9020 FORMAT(1X,A,F10.5,A)
668 9030 FORMAT(1X,A,F8.2,A)
669
670   END SUBROUTINE dia_fwb
671
672#else
673   !!----------------------------------------------------------------------
674   !!   Default option :                                       Dummy Module
675   !!----------------------------------------------------------------------
676   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .FALSE.    !: fresh water budget flag
677CONTAINS
678   SUBROUTINE dia_fwb( kt )        ! Empty routine
679      WRITE(*,*) 'dia_fwb: : You should not have seen this print! error?', kt
680   END SUBROUTINE dia_fwb
681#endif
682
683   !!======================================================================
684END MODULE diafwb