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.
diaar5.F90 in branches/UKMO/dev_r5518_GO6_fix_diaar5/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/dev_r5518_GO6_fix_diaar5/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90 @ 10870

Last change on this file since 10870 was 10870, checked in by frrh, 5 years ago

Apply correction to tnpeo calculation including reordering of k loop
as per #1884. This version of the code reflects the vn 3.6 STABLE branch
more than the vn 4.0 trunk. Although functionality is the same, here we
include the lbc_lnk of the STABLE branch (not used in the trunk)
and the workspace allocation of zpe as a pointer rather than the trunk
version which defines zpe (and other workspace arrays) simply as normal
allocatable arrays. i.e. I've opted to reflect the STABLE branch exactly
rather than to transplant more optimal code from the trunk out of context.
File size: 13.9 KB
Line 
1MODULE diaar5
2   !!======================================================================
3   !!                       ***  MODULE  diaar5  ***
4   !! AR5 diagnostics
5   !!======================================================================
6   !! History :  3.2  !  2009-11  (S. Masson)  Original code
7   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
8   !!----------------------------------------------------------------------
9#if defined key_diaar5   || defined key_esopa
10   !!----------------------------------------------------------------------
11   !!   'key_diaar5'  :                           activate ar5 diagnotics
12   !!----------------------------------------------------------------------
13   !!   dia_ar5       : AR5 diagnostics
14   !!   dia_ar5_init  : initialisation of AR5 diagnostics
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and active tracers
17   USE dom_oce        ! ocean space and time domain
18   USE eosbn2         ! equation of state                (eos_bn2 routine)
19   USE lib_mpp        ! distribued memory computing library
20   USE iom            ! I/O manager library
21   USE timing         ! preformance summary
22   USE wrk_nemo       ! working arrays
23   USE fldread        ! type FLD_N
24   USE phycst         ! physical constant
25   USE in_out_manager  ! I/O manager
26   USE zdfddm
27   USE zdf_oce
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   dia_ar5        ! routine called in step.F90 module
33   PUBLIC   dia_ar5_init   ! routine called in opa.F90 module
34   PUBLIC   dia_ar5_alloc  ! routine called in nemogcm.F90 module
35
36   LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .TRUE.   ! coupled flag
37
38   REAL(wp)                         ::   vol0         ! ocean volume (interior domain)
39   REAL(wp)                         ::   area_tot     ! total ocean surface (interior domain)
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   area         ! cell surface (interior domain)
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   thick0       ! ocean thickness (interior domain)
42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity
43     
44   !! * Substitutions
45#  include "domzgr_substitute.h90"
46#  include "zdfddm_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   FUNCTION dia_ar5_alloc()
55      !!----------------------------------------------------------------------
56      !!                    ***  ROUTINE dia_ar5_alloc  ***
57      !!----------------------------------------------------------------------
58      INTEGER :: dia_ar5_alloc
59      !!----------------------------------------------------------------------
60      !
61      ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc )
62      !
63      IF( lk_mpp             )   CALL mpp_sum ( dia_ar5_alloc )
64      IF( dia_ar5_alloc /= 0 )   CALL ctl_warn('dia_ar5_alloc: failed to allocate arrays')
65      !
66   END FUNCTION dia_ar5_alloc
67
68
69   SUBROUTINE dia_ar5( kt )
70      !!----------------------------------------------------------------------
71      !!                    ***  ROUTINE dia_ar5  ***
72      !!
73      !! ** Purpose :   compute and output some AR5 diagnostics
74      !!----------------------------------------------------------------------
75      !
76      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
77      !
78      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
79      REAL(wp) ::   zvolssh, zvol, zssh_steric, zztmp, zarho, ztemp, zsal, zmass
80      REAL(wp) ::   zaw, zbw, zrw
81      !
82      REAL(wp), POINTER, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace
83      REAL(wp), POINTER, DIMENSION(:,:) :: zpe                        ! 2D workspace
84      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zrhd , zrhop               ! 3D workspace
85      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace
86      !!--------------------------------------------------------------------
87      IF( nn_timing == 1 )   CALL timing_start('dia_ar5')
88
89      !Call to init moved to here so that we can call iom_use in the
90      !initialisation
91      IF( kt == nit000 )     CALL dia_ar5_init
92 
93      CALL wrk_alloc( jpi , jpj              , zarea_ssh , zbotpres, zpe)
94      CALL wrk_alloc( jpi , jpj , jpk        , zrhd      , zrhop    )
95      CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 )
96
97      zarea_ssh(:,:) = area(:,:) * sshn(:,:)
98
99      !                                         ! total volume of liquid seawater
100      zvolssh = SUM( zarea_ssh(:,:) ) 
101      IF( lk_mpp )   CALL mpp_sum( zvolssh )
102      zvol = vol0 + zvolssh
103     
104      CALL iom_put( 'voltot', zvol               )
105      CALL iom_put( 'sshtot', zvolssh / area_tot )
106      CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) )
107
108      !                     
109      IF( iom_use('sshthster')) THEN
110         ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh
111         ztsn(:,:,:,jp_sal) = sn0(:,:,:)
112         CALL eos( ztsn, zrhd, fsdept_n(:,:,:) )                       ! now in situ density using initial salinity
113         !
114         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
115         DO jk = 1, jpkm1
116            zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)
117         END DO
118         IF( .NOT.lk_vvl ) THEN
119            IF ( ln_isfcav ) THEN
120               DO ji=1,jpi
121                  DO jj=1,jpj
122                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
123                  END DO
124               END DO
125            ELSE
126               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
127            END IF
128         END IF
129      !                                         
130         zarho = SUM( area(:,:) * zbotpres(:,:) ) 
131         IF( lk_mpp )   CALL mpp_sum( zarho )
132         zssh_steric = - zarho / area_tot
133         CALL iom_put( 'sshthster', zssh_steric )
134      ENDIF
135     
136      !                                         ! steric sea surface height
137      CALL eos( tsn, zrhd, zrhop, fsdept_n(:,:,:) )                 ! now in situ and potential density
138      zrhop(:,:,jpk) = 0._wp
139      CALL iom_put( 'rhop', zrhop )
140      !
141      zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
142      DO jk = 1, jpkm1
143         zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)
144      END DO
145      IF( .NOT.lk_vvl ) THEN
146         IF ( ln_isfcav ) THEN
147            DO ji=1,jpi
148               DO jj=1,jpj
149                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
150               END DO
151            END DO
152         ELSE
153            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
154         END IF
155      END IF
156      !   
157      zarho = SUM( area(:,:) * zbotpres(:,:) ) 
158      IF( lk_mpp )   CALL mpp_sum( zarho )
159      zssh_steric = - zarho / area_tot
160      CALL iom_put( 'sshsteric', zssh_steric )
161     
162      !                                         ! ocean bottom pressure
163      zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
164      zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
165      CALL iom_put( 'botpres', zbotpres )
166
167      !                                         ! Mean density anomalie, temperature and salinity
168      ztemp = 0._wp
169      zsal  = 0._wp
170      DO jk = 1, jpkm1
171         DO jj = 1, jpj
172            DO ji = 1, jpi
173               zztmp = area(ji,jj) * fse3t(ji,jj,jk)
174               ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem)
175               zsal  = zsal  + zztmp * tsn(ji,jj,jk,jp_sal)
176            END DO
177         END DO
178      END DO
179      IF( .NOT.lk_vvl ) THEN
180         IF ( ln_isfcav ) THEN
181            DO ji=1,jpi
182               DO jj=1,jpj
183                  ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem) 
184                  zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal) 
185               END DO
186            END DO
187         ELSE
188            ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
189            zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
190         END IF
191      ENDIF
192      IF( lk_mpp ) THEN 
193         CALL mpp_sum( ztemp )
194         CALL mpp_sum( zsal  )
195      END IF
196      !
197      zmass = rau0 * ( zarho + zvol )                 ! total mass of liquid seawater
198      ztemp = ztemp / zvol                            ! potential temperature in liquid seawater
199      zsal  = zsal  / zvol                            ! Salinity of liquid seawater
200      !
201      CALL iom_put( 'masstot', zmass )
202      CALL iom_put( 'temptot', ztemp )
203      CALL iom_put( 'saltot' , zsal  )
204
205      IF( iom_use( 'tnpeo' )) THEN   
206      ! Work done against stratification by vertical mixing
207      ! Exclude points where rn2 is negative as convection kicks in here and
208      ! work is not being done against stratification
209          zpe(:,:) = 0._wp
210          IF( ln_zdfddm ) THEN
211             DO jk = 2, jpk
212                DO jj = 1, jpj
213                   DO ji = 1, jpi
214                     IF( rn2(ji,jj,jk) > 0._wp ) THEN
215
216                        zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   &
217                           &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) )
218                        !
219                        zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
220                        zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
221                        !
222
223                        zpe(ji, jj) = zpe(ji, jj)            &
224                                    -  grav * (  avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  &
225                                               - avs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )
226
227                   ENDDO
228                ENDDO
229             ENDDO
230          ELSE
231             DO jk=1,jpk
232                DO ji=1,jpi
233                   DO jj=1,jpj
234                      zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
235                   ENDDO
236                ENDDO
237             ENDDO
238          ENDIF
239          CALL lbc_lnk(zpe, 'T', 1._wp) 
240          CALL iom_put( 'tnpeo', zpe )         
241      ENDIF
242      !
243      CALL wrk_dealloc( jpi , jpj              , zarea_ssh , zbotpres, zpe )
244      CALL wrk_dealloc( jpi , jpj , jpk        , zrhd      , zrhop    )
245      CALL wrk_dealloc( jpi , jpj , jpk , jpts , ztsn                 )
246      !
247      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5')
248      !
249   END SUBROUTINE dia_ar5
250
251
252   SUBROUTINE dia_ar5_init
253      !!----------------------------------------------------------------------
254      !!                  ***  ROUTINE dia_ar5_init  ***
255      !!                   
256      !! ** Purpose :   initialization for AR5 diagnostic computation
257      !!----------------------------------------------------------------------
258      INTEGER  ::   inum
259      INTEGER  ::   ik
260      INTEGER  ::   ji, jj, jk  ! dummy loop indices
261      REAL(wp) ::   zztmp 
262      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity
263      !
264      !!----------------------------------------------------------------------
265      !
266      IF( nn_timing == 1 )   CALL timing_start('dia_ar5_init')
267      !
268      CALL wrk_alloc( jpi, jpj, jpk, 2, zsaldta )
269      !                                      ! allocate dia_ar5 arrays
270      IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
271
272      area(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
273
274      area_tot = SUM( area(:,:) )   ;   IF( lk_mpp )   CALL mpp_sum( area_tot )
275
276      vol0        = 0._wp
277      thick0(:,:) = 0._wp
278      DO jk = 1, jpkm1
279         vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) )
280         thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk)
281      END DO
282      IF( lk_mpp )   CALL mpp_sum( vol0 )
283
284      IF( iom_use('sshthster')) THEN
285         CALL iom_open ( 'sali_ref_clim_monthly', inum )
286         CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  )
287         CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 )
288         CALL iom_close( inum )
289
290         sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )       
291         sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
292         IF( ln_zps ) THEN               ! z-coord. partial steps
293            DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
294               DO ji = 1, jpi
295                  ik = mbkt(ji,jj)
296                  IF( ik > 1 ) THEN
297                     zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
298                     sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
299                  ENDIF
300               END DO
301            END DO
302         ENDIF
303      ENDIF
304      !
305      CALL wrk_dealloc( jpi, jpj, jpk, 2, zsaldta )
306      !
307      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init')
308      !
309   END SUBROUTINE dia_ar5_init
310
311#else
312   !!----------------------------------------------------------------------
313   !!   Default option :                                         NO diaar5
314   !!----------------------------------------------------------------------
315   LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .FALSE.   ! coupled flag
316CONTAINS
317   SUBROUTINE dia_ar5_init    ! Dummy routine
318   END SUBROUTINE dia_ar5_init
319   SUBROUTINE dia_ar5( kt )   ! Empty routine
320      INTEGER ::   kt
321      WRITE(*,*) 'dia_ar5: You should not have seen this print! error?', kt
322   END SUBROUTINE dia_ar5
323#endif
324
325   !!======================================================================
326END MODULE diaar5
Note: See TracBrowser for help on using the repository browser.