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

Ticket #127: diafwb_gm_style.F90

File diafwb_gm_style.F90, 18.0 KB (added by gm, 15 years ago)

bug correction for mpp from Gurvan + style and minor optimisation

Line 
1MODULE diafwb
2   !!======================================================================
3   !!                       ***  MODULE  diafwb  ***
4   !! Ocean diagnostics: freshwater budget + some straits
5   !!
6   !!                 ----   ORCA configuration only   ----
7   !!
8   !!======================================================================
9   !! History :  OPA  !  2001-02  (E. Durand)  Original code
10   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
11   !!            3.2  !  2009-07  (A.R. Porter, G. Madec) correct a mpp bug + salt content computation
12   !!----------------------------------------------------------------------
13#if ( defined key_orca_r2 || defined  key_orca_r4 ) && ! defined key_dynspg_rl && ! defined key_coupled
14   !!----------------------------------------------------------------------
15   !!   NOT "key_dynspg_rl" and "key_orca_r2 or 4"
16   !!----------------------------------------------------------------------
17   !!   dia_fwb     : freshwater budget for global ocean configurations
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and tracers
20   USE dom_oce         ! ocean space and time domain
21   USE phycst          ! physical constants
22   USE sbc_oce         ! ???
23   USE zdf_oce         ! ocean vertical physics
24   USE in_out_manager  ! I/O manager
25   USE lib_mpp         ! distributed memory computing library
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC dia_fwb    ! routine called by step.F90
31
32   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .TRUE.    !: fresh water budget flag
33
34   REAL(wp) ::   a_emp 
35   REAL(wp) ::   a_sshb, a_sshn, a_salb, a_saln
36   REAL(wp), DIMENSION(4) ::   a_flxi, a_temi, a_sali    ! ???
37   REAL(wp), DIMENSION(4) ::   a_flxo, a_temo, a_salo    ! ???
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
44   !! $Id$
45   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47
48CONTAINS
49
50   SUBROUTINE dia_fwb( kt )
51      !!---------------------------------------------------------------------
52      !!                  ***  ROUTINE dia_fwb  ***
53      !!     
54      !! ** Purpose :   ORCA diagnostics : salt content, volume and some straits
55      !!----------------------------------------------------------------------
56      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
57      !!
58      INTEGER  ::   inum             ! temporary logical unit
59      INTEGER  ::   ji, jj, jk, jt   ! dummy loop indices
60      INTEGER  ::   ii0, ii1, ij0, ij1
61      REAL(wp) ::   zarea, zvol, zwei
62      REAL(wp) ::   zt, zs, zu, zup, zum 
63      REAL(wp) ::   zsm0, zcoef
64      REAL(wp), DIMENSION(4)       ::  ztemi, ztemo, zsali, zsalo, zflxi, zflxo
65      REAL(wp), DIMENSION(jpi,jpj) ::   zsurf
66      !!----------------------------------------------------------------------
67
68      !                                ! ----------------------- !
69      IF( cp_cfg == "orca" ) THEN      !   ORCA configurations   !
70                                       ! ----------------------- !
71
72         ! Volume & Salt global mean diag.
73         ! ------------------------------
74         zsm0 = 34.72654      ! Mean global salinity
75
76         IF( kt == nit000 ) THEN      ! Initial values
77            !                                          ! interior masked surface
78            zsurf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:) 
79            !
80            a_sshb = SUM( zsurf(:,:) * sshb(:,:) )     ! mean initial ssh
81            IF( lk_mpp )   CALL mpp_sum( a_sshb )           ! sum over the global domain
82            !
83            a_salb   = 0.e0                            ! mean initial salinity
84            DO jk = 1, jpkm1
85               a_salb = a_salb + SUM(  zsurf(:,:) * fse3t(:,:,jk) * ( sb(:,:,jk) - zsm0 ) * tmask(:,:,jk)  )
86            END DO
87            IF( .NOT. lk_vvl )   a_salt = a_salt + SUM( zsurf(:,:) * ssh(:,:) * ( sb(:,:,1) - zsm0 ) )
88            IF( lk_mpp )   CALL mpp_sum( a_salb )           ! sum over the global domain
89            !
90            a_emp = 0.e0                               ! prepare the cumulation of emp
91         ENDIF
92     
93         a_emp = a_emp + SUM( zsurf(:,:) * emp(:,:) )  ! cumulative evap - precip - runoffs
94         IF( lk_mpp )   CALL mpp_sum( a_emp    )            ! sum over the global domain
95
96
97         IF( kt == nitend ) THEN      ! Final values
98            !                                          ! interior masked surface
99            zsurf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:) 
100            !
101            zarea  = SUM( zsurf(:,:) )                 ! surface area
102            IF( lk_mpp )   CALL mpp_sum( a_sshn )           ! sum over the global domain
103            !
104            a_sshn = SUM( zsurf(:,:) * sshn(:,:) )     ! Mean sea level at nitend
105            IF( lk_mpp )   CALL mpp_sum( zarea  )           ! sum over the global domain
106            !
107            a_saln = 0.e0                              ! Mean salinity and volume
108            zvol   = 0.e0
109            DO jk = 1, jpkm1   
110               DO jj = 2, jpjm1
111                  DO ji = fs_2, fs_jpim1   ! vector opt.
112                     zwei   = zsurf(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)
113                     a_saln = a_saln + ( sn(ji,jj,jk) - zsm0 ) * zwei
114                     zvol   = zvol  + zwei
115                  END DO
116               END DO
117            END DO
118            IF( .NOT. lk_vvl ) THEN                    ! correct with salt content in the ssh 'layer'
119               a_saln = a_saln + SUM(  zsurf(:,:) * ssh(:,:) * ( sb(:,:,1) - zsm0 )  )
120               zvol   = zvol   + SUM(  zsurf(:,:) * ssh(:,:)                         )
121            ENDIF
122            IF( lk_mpp )   CALL mpp_sum( a_saln )           ! sum over the global domain
123            IF( lk_mpp )   CALL mpp_sum( zvol   )           ! sum over the global domain
124            !
125            a_emp = a_emp * rdttra(1) * 1.e-3               ! Conversion in m3
126            !
127         ENDIF
128
129
130         ! Transport through some straits
131         ! ------------------------------
132         ! 1 --> Gibraltar             
133         ! 2 --> Cadiz
134         ! 3 --> Red Sea
135         ! 4 --> Baltic Sea
136         !
137         IF( kt == nit000 ) THEN
138            a_flxi(:) = 0.e0   ;   a_temi(:) = 0.e0   ;   a_sali(:) = 0.e0
139            a_flxo(:) = 0.e0   ;   a_temo(:) = 0.e0   ;   a_salo(:) = 0.e0
140         ENDIF
141         zflxi(:) = 0.e0   ;   ztemi(:) = 0.e0   ;   zsali(:) = 0.e0
142         zflxo(:) = 0.e0   ;   ztemo(:) = 0.e0   ;   zsalo(:) = 0.e0
143         !
144         !* Mean flow at Gibraltar (1)
145         SELECT CASE ( jp_cfg )                 ! select Gribraltar position
146         CASE ( 4 )                                  ! ORCA_R4
147            ii0 = 70   ;   ii1 = 70
148            ij0 = 52   ;   ij1 = 52
149         CASE ( 2 )                                  ! ORCA_R2
150            ii0 = 140   ;   ii1 = 140
151            ij0 = 102   ;   ij1 = 102
152         CASE DEFAULT                                ! other not yet implemented
153            CALL ctl_stop( ' dia_fwb Not yet implemented for ORCA_R1, _R05 or R025' )
154         END SELECT
155         !
156         DO jj = mj0(ij0), mj1(ij1)
157            DO ji = mi0(ii0), MAX( mi1(ii1), jpim1 )      ! max used to avoid mpp out-of-bound
158               DO jk = 1, jpkm1
159                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
160                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
161                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)   ! only if interior point
162                  zup = 0.5 * ( zu + ABS(zu) )
163                  zum = 0.5 * ( zu - ABS(zu) )
164                  !
165                  zflxi(1) = zflxi(1) +      zup
166                  ztemi(1) = ztemi(1) + zt * zup
167                  zsali(1) = zsali(1) + zs * zup
168                  !
169                  zflxo(1) = zflxo(1) +      zum
170                  ztemo(1) = ztemo(1) + zt * zum
171                  zsalo(1) = zsalo(1) + zs * zum
172               END DO
173            END DO
174         END DO
175         !     
176         !* Mean flow at Cadiz (2)
177         SELECT CASE ( jp_cfg )                 ! select Cadiz entrance position
178         CASE ( 4 )                                  ! ORCA_R4 configuration
179            ii0 = 69   ;   ii1 = 69
180            ij0 = 52   ;   ij1 = 52
181         CASE ( 2 )                                  ! ORCA_R2 configuration
182            ii0 = 137   ;   ii1 = 137
183            ij0 = 101   ;   ij1 = 102
184         CASE DEFAULT                                ! other not yet implemented
185            CALL ctl_stop( ' dia_fwb Not yet implemented for ORCA_R1, _R05 or R025' )
186         END SELECT
187         !
188         DO jj = mj0(ij0), mj1(ij1)
189            DO ji = mi0(ii0), MAX( mi1(ii1), jpim1 )      ! max used to avoid mpp out-of-bound
190               DO jk = 1, jpkm1 
191                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
192                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
193                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)   ! only if interior point
194                  zup = 0.5 * ( zu + ABS(zu) )
195                  zum = 0.5 * ( zu - ABS(zu) )
196                  !
197                  zflxi(2) = zflxi(2) +      zup
198                  ztemi(2) = ztemi(2) + zt * zup
199                  zsali(2) = zsali(2) + zs * zup
200                  !
201                  zflxo(2) = zflxo(2) +      zum
202                  ztemo(2) = ztemo(2) + zt * zum
203                  zsalo(2) = zsalo(2) + zs * zum
204               END DO
205            END DO
206         END DO
207
208         !* Bab-el-Mandeb (Red Sea entrance) (3)             
209         SELECT CASE ( jp_cfg )
210         CASE ( 4 )                                  !  ORCA_R4 configuration
211            ii0 = 83   ;   ii1 = 83
212            ij0 = 45   ;   ij1 = 45
213         CASE ( 2 )                                  !  ORCA_R2 configuration
214            ii0 = 160   ;   ii1 = 160
215            ij0 = 88    ;   ij1 = 88 
216         CASE DEFAULT                                ! other not yet implemented
217            CALL ctl_stop( ' dia_fwb Not yet implemented for ORCA_R1, _R05 or R025' )
218         END SELECT
219         !
220         DO jj = mj0(ij0), mj1(ij1)
221            DO ji = mi0(ii0), MAX( mi1(ii1), jpim1 )      ! max used to avoid mpp out-of-bound
222               DO jk = 1, jpk 
223                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
224                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
225                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)   ! only if interior point
226                  zup = 0.5 * ( zu + ABS(zu) )
227                  zum = 0.5 * ( zu - ABS(zu) )
228                  !
229                  zflxi(3) = zflxi(3) +      zup
230                  ztemi(3) = ztemi(3) + zt * zup
231                  zsali(3) = zsali(3) + zs * zup
232                  !
233                  zflxo(3) = zflxo(3) +      zum
234                  ztemo(3) = ztemo(3) + zt * zum
235                  zsalo(3) = zsalo(3) + zs * zum
236               END DO
237            END DO
238         END DO
239
240         !* Mean flow at Baltic Sea entrance (4)               
241         SELECT CASE ( jp_cfg )
242         CASE ( 4 )                                  !  ORCA_R4 configuration
243            ii0 = 1     ;   ii1 = 1 
244            ij0 = 1     ;   ij1 = 1 
245         CASE ( 2 )                                  !  ORCA_R2 configuration
246            ii0 = 146   ;   ii1 = 146 
247            ij0 = 116   ;   ij1 = 116
248         CASE DEFAULT                                ! other not yet implemented
249            CALL ctl_stop( ' dia_fwb Not yet implemented for ORCA_R1, _R05 or R025' )
250         END SELECT
251         !
252         DO jj = mj0(ij0), mj1(ij1)
253            DO ji = mi0(ii0), MAX( mi1(ii1), jpim1 )      ! max used to avoid mpp out-of-bound
254               DO jk = 1, jpk
255                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
256                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
257                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)   ! only if interior point
258                  zup = 0.5 * ( zu + ABS(zu) )
259                  zum = 0.5 * ( zu - ABS(zu) )
260                  !
261                  zflxi(4) = zflxi(4) +      zup
262                  ztemi(4) = ztemi(4) + zt * zup
263                  zsali(4) = zsali(4) + zs * zup
264                  !
265                  zflxo(4) = zflxo(4) +      zum
266                  ztemo(4) = ztemo(4) + zt * zum
267                  zsalo(4) = zsalo(4) + zs * zum
268               END DO
269            END DO
270         END DO
271         !
272         ! Time cumulation
273         DO jt = 1, 4 
274            IF( zflxi(jt) /= 0.e0 ) THEN
275               a_flxi(jt) = a_flxi(jt) + zflxi(jt)
276               a_temi(jt) = a_temi(jt) + ztemi(jt) / zflxi(jt)
277               a_sali(jt) = a_sali(jt) + zsali(jt) / zflxi(jt)
278            ENDIF
279            IF( zflxo(jt) /= 0.e0 ) THEN
280               a_flxo(jt) = a_flxo(jt) + zflxo(jt)
281               a_temo(jt) = a_temo(jt) + ztemo(jt) / zflxo(jt)
282               a_salo(jt) = a_salo(jt) + zsalo(jt) / zflxo(jt)
283            ENDIF
284         END DO
285
286         IF( kt == nitend ) THEN              ! time averaged
287            zcoef = 1.e0 / REAL( nitend - nit000 + 1 )
288            DO jt = 1, 4 
289               a_flxi(jt) = a_flxi(jt) * zcoef * 1.e-6   ;   a_flxo(jt) = a_flxo(jt) * zcoef * 1.e-6
290               a_temi(jt) = a_temi(jt) * zcoef           ;   a_temo(jt) = a_temo(jt) * zcoef
291               a_sali(jt) = a_sali(jt) * zcoef           ;   a_salo(jt) = a_salo(jt) * zcoef
292            END DO
293            IF( lk_mpp ) THEN                 ! sum over the global domain
294               CALL mpp_sum( a_flxi, 4 )                 ;   CALL mpp_sum( a_flxo, 4 ) 
295               CALL mpp_sum( a_temi, 4 )                 ;   CALL mpp_sum( a_temo, 4 )
296               CALL mpp_sum( a_sali, 4 )                 ;   CALL mpp_sum( a_salo, 4 )
297            ENDIF
298         ENDIF
299
300
301         ! Write the diagnostics
302         ! --------------------------
303         !
304         IF ( kt == nitend ) THEN
305            !
306            CALL ctlopn( inum, 'STRAIT.dat', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 )
307            WRITE(inum,*)
308            WRITE(inum,*)    'Net freshwater budget '
309            WRITE(inum,9010) '  emp    = ',a_emp,   ' m3 =', a_emp / (REAL(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
310            WRITE(inum,*)
311            WRITE(inum,9010) '  zarea =',zarea
312            WRITE(inum,9010) '  zvol  =',zvol
313            WRITE(inum,*)
314            WRITE(inum,*)    'Mean sea level : '
315            WRITE(inum,9010) '  at nit000 = ', a_sshb           , ' m3 '
316            WRITE(inum,9010) '  at nitend = ', a_sshn           , ' m3 '
317            WRITE(inum,9010) '  diff      = ', (a_sshn - a_sshb), ' m3 =', (a_sshn-a_sshb)/(REAL(nitend-nit000+1)*rdt) * 1.e-6,' Sv'
318            WRITE(inum,9020) '  mean sea level elevation    ='  , a_sshn / zarea,' m'
319            WRITE(inum,*)
320            WRITE(inum,*)    'Anomaly of salinity content : '
321            WRITE(inum,9010) '  at nit000 = ',  a_salb          , ' psu.m3 '
322            WRITE(inum,9010) '  at nitend = ',  a_saln          , ' psu.m3 '
323            WRITE(inum,9010) '  diff      = ', (a_saln - a_salb), ' psu.m3'
324            WRITE(inum,*)
325            WRITE(inum,*)    'Mean salinity : '
326            WRITE(inum,9020) '  at nit000 =', a_salb / zvol + zsm0    , ' psu '
327            WRITE(inum,9020) '  at nitend =', a_saln / zvol + zsm0    , ' psu '
328            WRITE(inum,9020) '  diff      =', (a_saln - a_salb) / zvol, ' psu'
329            WRITE(inum,9020) '  S-SLevitus=', a_saln / zvol           ,' psu'
330            WRITE(inum,*)
331            WRITE(inum,*)    'Gibraltar : '
332            WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(1)
333            WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(1)
334            WRITE(inum,9030) '  T entrant (deg)   :', a_temi(1)
335            WRITE(inum,9030) '  T sortant (deg)   :', a_temo(1)
336            WRITE(inum,9030) '  S entrant (psu)   :', a_sali(1)
337            WRITE(inum,9030) '  S sortant (psu)   :', a_salo(1)
338            WRITE(inum,*)
339            WRITE(inum,*)    'Cadiz : '
340            WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(2)
341            WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(2)
342            WRITE(inum,9030) '  T entrant (deg)   :', a_temi(2)
343            WRITE(inum,9030) '  T sortant (deg)   :', a_temo(2)
344            WRITE(inum,9030) '  S entrant (psu)   :', a_sali(2)
345            WRITE(inum,9030) '  S sortant (psu)   :', a_salo(2)
346            WRITE(inum,*)
347            WRITE(inum,*)    'Bab el Mandeb : '
348            WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(3)
349            WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(3)
350            WRITE(inum,9030) '  T entrant (deg)   :', a_temi(3)
351            WRITE(inum,9030) '  T sortant (deg)   :', a_temo(3)
352            WRITE(inum,9030) '  S entrant (psu)   :', a_sali(3)
353            WRITE(inum,9030) '  S sortant (psu)   :', a_salo(3)
354            WRITE(inum,*)
355            WRITE(inum,*)    'Baltic : '
356            WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(4)
357            WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(4)
358            WRITE(inum,9030) '  T entrant (deg)   :', a_temi(4)
359            WRITE(inum,9030) '  T sortant (deg)   :', a_temo(4)
360            WRITE(inum,9030) '  S entrant (psu)   :', a_sali(4) 
361            WRITE(inum,9030) '  S sortant (psu)   :', a_salo(4)
362            CLOSE(inum)
363         ENDIF
364         !
365      ENDIF
366
367 9005 FORMAT(1X,A,ES24.16)
368 9010 FORMAT(1X,A,ES12.5,A,F10.5,A)
369 9020 FORMAT(1X,A,F10.5,A)
370 9030 FORMAT(1X,A,F9.4,A)
371      !
372   END SUBROUTINE dia_fwb
373
374#else
375   !!----------------------------------------------------------------------
376   !!   Default option :                                       Dummy Module
377   !!----------------------------------------------------------------------
378   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .FALSE.    !: fresh water budget flag
379CONTAINS
380   SUBROUTINE dia_fwb( kt )        ! Empty routine
381      WRITE(*,*) 'dia_fwb: : You should not have seen this print! error?', kt
382   END SUBROUTINE dia_fwb
383#endif
384
385   !!======================================================================
386END MODULE diafwb