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

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 19.8 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
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            DO jj = mj0(ij0), mj1(ij1)
231               DO jk = 1, jpk 
232                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
233                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
234                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
235                 
236                  IF( un(ji,jj,jk) > 0.e0 ) THEN
237                     zflxi(2) = zflxi(2) +    zu
238                     ztemi(2) = ztemi(2) + zt*zu
239                     zsali(2) = zsali(2) + zs*zu
240                  ELSE
241                     zflxo(2) = zflxo(2) +    zu
242                     ztemo(2) = ztemo(2) + zt*zu
243                     zsalo(2) = zsalo(2) + zs*zu
244                  ENDIF
245               END DO
246            END DO
247         END DO
248      ENDIF
249
250      ! Mean flow at Red Sea entrance
251      IF( cp_cfg == "orca" ) THEN
252               
253         SELECT CASE ( jp_cfg )
254         !                                           ! =======================
255         CASE ( 4 )                                  !  ORCA_R4 configuration
256            !                                        ! =======================
257            ii0 = 83   ;   ii1 = 83
258            ij0 = 45   ;   ij1 = 45
259            !                                        ! =======================
260         CASE ( 2 )                                  !  ORCA_R2 configuration
261            !                                        ! =======================
262            ii0 = 160   ;   ii1 = 160
263            ij0 = 88    ;   ij1 = 88 
264            !                                        ! =======================
265         CASE ( 1 )                                  !  ORCA_R1 configurations
266            !                                        ! =======================
267            ii0 = 331   ;   ii1 = 331
268            ij0 = 176   ;   ij1 = 176
269            !                                        ! =======================
270         CASE DEFAULT                                !    ORCA R05 or R025
271            !                                        ! =======================
272            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
273            !
274         END SELECT
275         !
276         DO ji = mi0(ii0), mi1(ii1)
277            DO jj = mj0(ij0), mj1(ij1)
278               DO jk = 1, jpk 
279                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
280                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
281                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
282                 
283                  IF( un(ji,jj,jk) > 0.e0 ) THEN
284                     zflxi(3) = zflxi(3) +    zu
285                     ztemi(3) = ztemi(3) + zt*zu
286                     zsali(3) = zsali(3) + zs*zu
287                  ELSE
288                     zflxo(3) = zflxo(3) +    zu
289                     ztemo(3) = ztemo(3) + zt*zu
290                     zsalo(3) = zsalo(3) + zs*zu
291                  ENDIF
292               END DO
293            END DO
294         END DO
295      ENDIF
296
297      ! Mean flow at Baltic Sea entrance
298      IF( cp_cfg == "orca" ) THEN
299               
300         SELECT CASE ( jp_cfg )
301         !                                           ! =======================
302         CASE ( 4 )                                  !  ORCA_R4 configuration
303            !                                        ! =======================
304            ii0 = 1     ;   ii1 = 1 
305            ij0 = 1     ;   ij1 = 1 
306            !                                        ! =======================
307         CASE ( 2 )                                  !  ORCA_R2 configuration
308            !                                        ! =======================
309            ii0 = 146   ;   ii1 = 146 
310            ij0 = 116   ;   ij1 = 116
311            !                                        ! =======================
312         CASE ( 1 )                                  !  ORCA_R1 configurations
313            !                                        ! =======================
314            ii0 = 297   ;   ii1 = 297 
315            ij0 = 230   ;   ij1 = 230
316            !                                        ! =======================
317         CASE DEFAULT                                !    ORCA R05 or R025
318            !                                        ! =======================
319            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
320            !
321         END SELECT
322         !
323         DO ji = mi0(ii0), mi1(ii1)
324            DO jj = mj0(ij0), mj1(ij1)
325               DO jk = 1, jpk
326                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
327                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
328                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
329                 
330                  IF( un(ji,jj,jk) > 0.e0 ) THEN
331                     zflxi(4) = zflxi(4) +    zu
332                     ztemi(4) = ztemi(4) + zt*zu
333                     zsali(4) = zsali(4) + zs*zu
334                  ELSE
335                     zflxo(4) = zflxo(4) +    zu
336                     ztemo(4) = ztemo(4) + zt*zu
337                     zsalo(4) = zsalo(4) + zs*zu
338                  ENDIF
339               END DO
340            END DO
341         END DO
342      ENDIF
343
344      ! Sum at each time-step
345      DO jt = 1, 4 
346         !
347         IF( zflxi(jt) /= 0.e0 ) THEN
348            a_flxi(jt) = a_flxi(jt) + zflxi(jt)
349            a_temi(jt) = a_temi(jt) + ztemi(jt)/zflxi(jt)
350            a_sali(jt) = a_sali(jt) + zsali(jt)/zflxi(jt)
351         ENDIF
352         !
353         IF( zflxo(jt) /= 0.e0 ) THEN
354            a_flxo(jt) = a_flxo(jt) + zflxo(jt)
355            a_temo(jt) = a_temo(jt) + ztemo(jt)/zflxo(jt)
356            a_salo(jt) = a_salo(jt) + zsalo(jt)/zflxo(jt)
357         ENDIF
358         !
359      END DO
360
361      IF( kt == nitend ) THEN
362         DO jt = 1, 4 
363            a_flxi(jt) = a_flxi(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
364            a_temi(jt) = a_temi(jt) /   FLOAT( nitend - nit000 + 1 )
365            a_sali(jt) = a_sali(jt) /   FLOAT( nitend - nit000 + 1 )
366            a_flxo(jt) = a_flxo(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
367            a_temo(jt) = a_temo(jt) /   FLOAT( nitend - nit000 + 1 )
368            a_salo(jt) = a_salo(jt) /   FLOAT( nitend - nit000 + 1 )
369         END DO
370         IF( lk_mpp ) THEN
371            CALL mpp_sum( a_flxi, 4 )      ! sum over the global domain
372            CALL mpp_sum( a_temi, 4 )      ! sum over the global domain
373            CALL mpp_sum( a_sali, 4 )      ! sum over the global domain
374
375            CALL mpp_sum( a_flxo, 4 )      ! sum over the global domain
376            CALL mpp_sum( a_temo, 4 )      ! sum over the global domain
377            CALL mpp_sum( a_salo, 4 )      ! sum over the global domain
378         ENDIF
379      ENDIF
380
381
382      ! Ecriture des diagnostiques
383      ! --------------------------
384
385      IF ( kt == nitend .AND. cp_cfg == "orca" .AND. lwp ) THEN
386
387         CALL ctl_opn( inum, 'STRAIT.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
388         WRITE(inum,*)
389         WRITE(inum,*)    'Net freshwater budget '
390         WRITE(inum,9010) '  fwf    = ',a_fwf,   ' m3 =', a_fwf   /(FLOAT(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
391         WRITE(inum,*)
392         WRITE(inum,9010) '  zarea =',zarea
393         WRITE(inum,9010) '  zvol  =',zvol
394         WRITE(inum,*)
395         WRITE(inum,*)    'Mean sea level : '
396         WRITE(inum,9010) '  at nit000 = ',a_sshb        ,' m3 '
397         WRITE(inum,9010) '  at nitend = ',a_sshn        ,' m3 '
398         WRITE(inum,9010) '  diff      = ',(a_sshn-a_sshb),' m3 =', (a_sshn-a_sshb)/(FLOAT(nitend-nit000+1)*rdt) * 1.e-6,' Sv'
399         WRITE(inum,9020) '  mean sea level elevation    =', a_sshn/zarea,' m'
400         WRITE(inum,*)
401         WRITE(inum,*)    'Anomaly of salinity content : '
402         WRITE(inum,9010) '  at nit000 = ',a_salb        ,' psu.m3 '
403         WRITE(inum,9010) '  at nitend = ',a_saln        ,' psu.m3 '
404         WRITE(inum,9010) '  diff      = ',(a_saln-a_salb),' psu.m3'
405         WRITE(inum,*)
406         WRITE(inum,*)    'Mean salinity : '
407         WRITE(inum,9020) '  at nit000 =',a_salb/zvol+zsm0   ,' psu '
408         WRITE(inum,9020) '  at nitend =',a_saln/zvol+zsm0   ,' psu '
409         WRITE(inum,9020) '  diff      =',(a_saln-a_salb)/zvol,' psu'
410         WRITE(inum,9020) '  S-SLevitus=',a_saln/zvol,' psu'
411         WRITE(inum,*)
412         WRITE(inum,*)    'Gibraltar : '
413         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(1)
414         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(1)
415         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(1)
416         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(1)
417         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(1)
418         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(1)
419         WRITE(inum,*)
420         WRITE(inum,*)    'Cadiz : '
421         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(2)
422         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(2)
423         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(2)
424         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(2)
425         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(2)
426         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(2)
427         WRITE(inum,*)
428         WRITE(inum,*)    'Bab el Mandeb : '
429         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(3)
430         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(3)
431         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(3)
432         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(3)
433         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(3)
434         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(3)
435         WRITE(inum,*)
436         WRITE(inum,*)    'Baltic : '
437         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(4)
438         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(4)
439         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(4)
440         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(4)
441         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(4)
442         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(4)
443         CLOSE(inum)
444      ENDIF
445
446 9005 FORMAT(1X,A,ES24.16)
447 9010 FORMAT(1X,A,ES12.5,A,F10.5,A)
448 9020 FORMAT(1X,A,F10.5,A)
449 9030 FORMAT(1X,A,F9.4,A)
450
451   END SUBROUTINE dia_fwb
452
453#else
454   !!----------------------------------------------------------------------
455   !!   Default option :                                       Dummy Module
456   !!----------------------------------------------------------------------
457   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .FALSE.    !: fresh water budget flag
458CONTAINS
459   SUBROUTINE dia_fwb( kt )        ! Empty routine
460      WRITE(*,*) 'dia_fwb: : You should not have seen this print! error?', kt
461   END SUBROUTINE dia_fwb
462#endif
463
464   !!======================================================================
465END MODULE diafwb
Note: See TracBrowser for help on using the repository browser.