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/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90 @ 5075

Last change on this file since 5075 was 5075, checked in by timgraham, 9 years ago

Upgraded branch to current head of trunk (r5072) so it can be used with the trunk

  • Property svn:keywords set to Id
File size: 19.5 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   !!----------------------------------------------------------------------
11   !!   Only for ORCA2 ORCA1 and ORCA025
12   !!----------------------------------------------------------------------
13   !!----------------------------------------------------------------------
14   !!   dia_fwb     : freshwater budget for global ocean configurations
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE phycst          ! physical constants
19   USE sbc_oce         ! ???
20   USE zdf_oce         ! ocean vertical physics
21   USE in_out_manager  ! I/O manager
22   USE lib_mpp         ! distributed memory computing library
23   USE timing          ! preformance summary
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC dia_fwb    ! routine called by step.F90
29
30   REAL(wp)               ::   a_fwf ,          &
31      &                        a_sshb, a_sshn, a_salb, a_saln
32   REAL(wp), DIMENSION(4) ::   a_flxi, a_flxo, a_temi, a_temo, a_sali, a_salo
33
34   !! * Substitutions
35#  include "domzgr_substitute.h90"
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
39   !! $Id$
40   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42
43CONTAINS
44
45   SUBROUTINE dia_fwb( kt )
46      !!---------------------------------------------------------------------
47      !!                  ***  ROUTINE dia_fwb  ***
48      !!     
49      !! ** Purpose :
50      !!----------------------------------------------------------------------
51      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
52      !!
53      INTEGER :: inum             ! temporary logical unit
54      INTEGER :: ji, jj, jk, jt   ! dummy loop indices
55      INTEGER :: ii0, ii1, ij0, ij1
56      REAL(wp) ::   zarea, zvol, zwei
57      REAL(wp) ::  ztemi(4), ztemo(4), zsali(4), zsalo(4), zflxi(4), zflxo(4)
58      REAL(wp) ::  zt, zs, zu 
59      REAL(wp) ::  zsm0, zfwfnew
60      IF( cp_cfg == "orca" .AND. jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) THEN
61      !!----------------------------------------------------------------------
62      IF( nn_timing == 1 )   CALL timing_start('dia_fwb')
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 + ( tsb(ji,jj,jk,jp_sal) - 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 + ( tsn(ji,jj,jk,jp_sal) - 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 ( 1 )                                  !  ORCA_R1 configurations
166            !                                        ! =======================
167            ii0 = 283   ;   ii1 = 283
168            ij0 = 200   ;   ij1 = 200
169            !                                        ! =======================
170         CASE DEFAULT                                !    ORCA R05 or R025
171            !                                        ! =======================
172            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
173            !
174         END SELECT
175         !
176         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
177            DO jj = mj0(ij0), mj1(ij1)
178               DO jk = 1, jpk 
179                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
180                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
181                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
182
183                  IF( un(ji,jj,jk) > 0.e0 ) THEN
184                     zflxi(1) = zflxi(1) +    zu
185                     ztemi(1) = ztemi(1) + zt*zu
186                     zsali(1) = zsali(1) + zs*zu
187                  ELSE
188                     zflxo(1) = zflxo(1) +    zu
189                     ztemo(1) = ztemo(1) + zt*zu
190                     zsalo(1) = zsalo(1) + zs*zu
191                  ENDIF
192               END DO
193            END DO
194         END DO
195      ENDIF
196     
197      ! Mean flow at Cadiz
198      IF( cp_cfg == "orca" ) THEN
199               
200         SELECT CASE ( jp_cfg )
201         !                                           ! =======================
202         CASE ( 4 )                                  !  ORCA_R4 configuration
203            !                                        ! =======================
204            ii0 = 69   ;   ii1 = 69
205            ij0 = 52   ;   ij1 = 52
206            !                                        ! =======================
207         CASE ( 2 )                                  !  ORCA_R2 configuration
208            !                                        ! =======================
209            ii0 = 137   ;   ii1 = 137
210            ij0 = 101   ;   ij1 = 102
211            !                                        ! =======================
212         CASE ( 1 )                                  !  ORCA_R1 configurations
213            !                                        ! =======================
214            ii0 = 282   ;   ii1 = 282
215            ij0 = 200   ;   ij1 = 200
216            !                                        ! =======================
217         CASE DEFAULT                                !    ORCA R05 or R025
218            !                                        ! =======================
219            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
220            !
221         END SELECT
222         !
223         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
224            DO jj = mj0(ij0), mj1(ij1)
225               DO jk = 1, jpk 
226                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
227                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
228                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
229                 
230                  IF( un(ji,jj,jk) > 0.e0 ) THEN
231                     zflxi(2) = zflxi(2) +    zu
232                     ztemi(2) = ztemi(2) + zt*zu
233                     zsali(2) = zsali(2) + zs*zu
234                  ELSE
235                     zflxo(2) = zflxo(2) +    zu
236                     ztemo(2) = ztemo(2) + zt*zu
237                     zsalo(2) = zsalo(2) + zs*zu
238                  ENDIF
239               END DO
240            END DO
241         END DO
242      ENDIF
243
244      ! Mean flow at Red Sea entrance
245      IF( cp_cfg == "orca" ) THEN
246               
247         SELECT CASE ( jp_cfg )
248         !                                           ! =======================
249         CASE ( 4 )                                  !  ORCA_R4 configuration
250            !                                        ! =======================
251            ii0 = 83   ;   ii1 = 83
252            ij0 = 45   ;   ij1 = 45
253            !                                        ! =======================
254         CASE ( 2 )                                  !  ORCA_R2 configuration
255            !                                        ! =======================
256            ii0 = 160   ;   ii1 = 160
257            ij0 = 88    ;   ij1 = 88 
258            !                                        ! =======================
259         CASE ( 1 )                                  !  ORCA_R1 configurations
260            !                                        ! =======================
261            ii0 = 331   ;   ii1 = 331
262            ij0 = 176   ;   ij1 = 176
263            !                                        ! =======================
264         CASE DEFAULT                                !    ORCA R05 or R025
265            !                                        ! =======================
266            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
267            !
268         END SELECT
269         !
270         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
271            DO jj = mj0(ij0), mj1(ij1)
272               DO jk = 1, jpk 
273                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
274                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
275                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
276                 
277                  IF( un(ji,jj,jk) > 0.e0 ) THEN
278                     zflxi(3) = zflxi(3) +    zu
279                     ztemi(3) = ztemi(3) + zt*zu
280                     zsali(3) = zsali(3) + zs*zu
281                  ELSE
282                     zflxo(3) = zflxo(3) +    zu
283                     ztemo(3) = ztemo(3) + zt*zu
284                     zsalo(3) = zsalo(3) + zs*zu
285                  ENDIF
286               END DO
287            END DO
288         END DO
289      ENDIF
290
291      ! Mean flow at Baltic Sea entrance
292      IF( cp_cfg == "orca" ) THEN
293               
294         SELECT CASE ( jp_cfg )
295         !                                           ! =======================
296         CASE ( 4 )                                  !  ORCA_R4 configuration
297            !                                        ! =======================
298            ii0 = 1     ;   ii1 = 1 
299            ij0 = 1     ;   ij1 = 1 
300            !                                        ! =======================
301         CASE ( 2 )                                  !  ORCA_R2 configuration
302            !                                        ! =======================
303            ii0 = 146   ;   ii1 = 146 
304            ij0 = 116   ;   ij1 = 116
305            !                                        ! =======================
306         CASE ( 1 )                                  !  ORCA_R1 configurations
307            !                                        ! =======================
308            ii0 = 297   ;   ii1 = 297 
309            ij0 = 230   ;   ij1 = 230
310            !                                        ! =======================
311         CASE DEFAULT                                !    ORCA R05 or R025
312            !                                        ! =======================
313            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
314            !
315         END SELECT
316         !
317         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
318            DO jj = mj0(ij0), mj1(ij1)
319               DO jk = 1, jpk
320                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
321                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
322                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
323                 
324                  IF( un(ji,jj,jk) > 0.e0 ) THEN
325                     zflxi(4) = zflxi(4) +    zu
326                     ztemi(4) = ztemi(4) + zt*zu
327                     zsali(4) = zsali(4) + zs*zu
328                  ELSE
329                     zflxo(4) = zflxo(4) +    zu
330                     ztemo(4) = ztemo(4) + zt*zu
331                     zsalo(4) = zsalo(4) + zs*zu
332                  ENDIF
333               END DO
334            END DO
335         END DO
336      ENDIF
337
338      ! Sum at each time-step
339      DO jt = 1, 4 
340         !
341         IF( zflxi(jt) /= 0.e0 ) THEN
342            a_flxi(jt) = a_flxi(jt) + zflxi(jt)
343            a_temi(jt) = a_temi(jt) + ztemi(jt)/zflxi(jt)
344            a_sali(jt) = a_sali(jt) + zsali(jt)/zflxi(jt)
345         ENDIF
346         !
347         IF( zflxo(jt) /= 0.e0 ) THEN
348            a_flxo(jt) = a_flxo(jt) + zflxo(jt)
349            a_temo(jt) = a_temo(jt) + ztemo(jt)/zflxo(jt)
350            a_salo(jt) = a_salo(jt) + zsalo(jt)/zflxo(jt)
351         ENDIF
352         !
353      END DO
354
355      IF( kt == nitend ) THEN
356         DO jt = 1, 4 
357            a_flxi(jt) = a_flxi(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
358            a_temi(jt) = a_temi(jt) /   FLOAT( nitend - nit000 + 1 )
359            a_sali(jt) = a_sali(jt) /   FLOAT( nitend - nit000 + 1 )
360            a_flxo(jt) = a_flxo(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
361            a_temo(jt) = a_temo(jt) /   FLOAT( nitend - nit000 + 1 )
362            a_salo(jt) = a_salo(jt) /   FLOAT( nitend - nit000 + 1 )
363         END DO
364         IF( lk_mpp ) THEN
365            CALL mpp_sum( a_flxi, 4 )      ! sum over the global domain
366            CALL mpp_sum( a_temi, 4 )      ! sum over the global domain
367            CALL mpp_sum( a_sali, 4 )      ! sum over the global domain
368
369            CALL mpp_sum( a_flxo, 4 )      ! sum over the global domain
370            CALL mpp_sum( a_temo, 4 )      ! sum over the global domain
371            CALL mpp_sum( a_salo, 4 )      ! sum over the global domain
372         ENDIF
373      ENDIF
374
375
376      ! Ecriture des diagnostiques
377      ! --------------------------
378
379      IF ( kt == nitend .AND. cp_cfg == "orca" .AND. lwp ) THEN
380
381         CALL ctl_opn( inum, 'STRAIT.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
382         WRITE(inum,*)
383         WRITE(inum,*)    'Net freshwater budget '
384         WRITE(inum,9010) '  fwf    = ',a_fwf,   ' m3 =', a_fwf   /(FLOAT(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
385         WRITE(inum,*)
386         WRITE(inum,9010) '  zarea =',zarea
387         WRITE(inum,9010) '  zvol  =',zvol
388         WRITE(inum,*)
389         WRITE(inum,*)    'Mean sea level : '
390         WRITE(inum,9010) '  at nit000 = ',a_sshb        ,' m3 '
391         WRITE(inum,9010) '  at nitend = ',a_sshn        ,' m3 '
392         WRITE(inum,9010) '  diff      = ',(a_sshn-a_sshb),' m3 =', (a_sshn-a_sshb)/(FLOAT(nitend-nit000+1)*rdt) * 1.e-6,' Sv'
393         WRITE(inum,9020) '  mean sea level elevation    =', a_sshn/zarea,' m'
394         WRITE(inum,*)
395         WRITE(inum,*)    'Anomaly of salinity content : '
396         WRITE(inum,9010) '  at nit000 = ',a_salb        ,' psu.m3 '
397         WRITE(inum,9010) '  at nitend = ',a_saln        ,' psu.m3 '
398         WRITE(inum,9010) '  diff      = ',(a_saln-a_salb),' psu.m3'
399         WRITE(inum,*)
400         WRITE(inum,*)    'Mean salinity : '
401         WRITE(inum,9020) '  at nit000 =',a_salb/zvol+zsm0   ,' psu '
402         WRITE(inum,9020) '  at nitend =',a_saln/zvol+zsm0   ,' psu '
403         WRITE(inum,9020) '  diff      =',(a_saln-a_salb)/zvol,' psu'
404         WRITE(inum,9020) '  S-SLevitus=',a_saln/zvol,' psu'
405         WRITE(inum,*)
406         WRITE(inum,*)    'Gibraltar : '
407         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(1)
408         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(1)
409         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(1)
410         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(1)
411         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(1)
412         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(1)
413         WRITE(inum,*)
414         WRITE(inum,*)    'Cadiz : '
415         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(2)
416         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(2)
417         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(2)
418         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(2)
419         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(2)
420         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(2)
421         WRITE(inum,*)
422         WRITE(inum,*)    'Bab el Mandeb : '
423         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(3)
424         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(3)
425         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(3)
426         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(3)
427         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(3)
428         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(3)
429         WRITE(inum,*)
430         WRITE(inum,*)    'Baltic : '
431         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(4)
432         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(4)
433         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(4)
434         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(4)
435         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(4)
436         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(4)
437         CLOSE(inum)
438      ENDIF
439
440      IF( nn_timing == 1 )   CALL timing_start('dia_fwb')
441
442 9005 FORMAT(1X,A,ES24.16)
443 9010 FORMAT(1X,A,ES12.5,A,F10.5,A)
444 9020 FORMAT(1X,A,F10.5,A)
445 9030 FORMAT(1X,A,F9.4,A)
446 
447      ENDIF
448
449   END SUBROUTINE dia_fwb
450
451   !!======================================================================
452END MODULE diafwb
Note: See TracBrowser for help on using the repository browser.