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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90 @ 4408

Last change on this file since 4408 was 4408, checked in by trackstand2, 10 years ago

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