source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90 @ 11035

Last change on this file since 11035 was 11035, checked in by frrh, 19 months ago

Apply changes from UKMO GMED ticket 453 to fix OOB bugs in
the calculation of diagnostic tnpeo and to apply the correct calculation
as per the current NEMO vn4.0 and 3.6_STABLE code bases, in place of the
existing incorrect formula. Thus, tnpeo diagnostic values will change
following this update.
Model prognostics and eveolution are unaffected.

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( lk_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                                               - fsavs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ))
226
227                      ENDIF
228                   ENDDO
229                ENDDO
230             ENDDO
231          ELSE
232             DO jk=1,jpk
233                DO ji=1,jpi
234                   DO jj=1,jpj
235                      zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
236                   ENDDO
237                ENDDO
238             ENDDO
239          ENDIF
240          CALL lbc_lnk(zpe, 'T', 1._wp) 
241          CALL iom_put( 'tnpeo', zpe )         
242      ENDIF
243      !
244      CALL wrk_dealloc( jpi , jpj              , zarea_ssh , zbotpres, zpe )
245      CALL wrk_dealloc( jpi , jpj , jpk        , zrhd      , zrhop    )
246      CALL wrk_dealloc( jpi , jpj , jpk , jpts , ztsn                 )
247      !
248      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5')
249      !
250   END SUBROUTINE dia_ar5
251
252
253   SUBROUTINE dia_ar5_init
254      !!----------------------------------------------------------------------
255      !!                  ***  ROUTINE dia_ar5_init  ***
256      !!                   
257      !! ** Purpose :   initialization for AR5 diagnostic computation
258      !!----------------------------------------------------------------------
259      INTEGER  ::   inum
260      INTEGER  ::   ik
261      INTEGER  ::   ji, jj, jk  ! dummy loop indices
262      REAL(wp) ::   zztmp 
263      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity
264      !
265      !!----------------------------------------------------------------------
266      !
267      IF( nn_timing == 1 )   CALL timing_start('dia_ar5_init')
268      !
269      CALL wrk_alloc( jpi, jpj, jpk, 2, zsaldta )
270      !                                      ! allocate dia_ar5 arrays
271      IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
272
273      area(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
274
275      area_tot = SUM( area(:,:) )   ;   IF( lk_mpp )   CALL mpp_sum( area_tot )
276
277      vol0        = 0._wp
278      thick0(:,:) = 0._wp
279      DO jk = 1, jpkm1
280         vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) )
281         thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk)
282      END DO
283      IF( lk_mpp )   CALL mpp_sum( vol0 )
284
285      IF( iom_use('sshthster')) THEN
286         CALL iom_open ( 'sali_ref_clim_monthly', inum )
287         CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  )
288         CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 )
289         CALL iom_close( inum )
290
291         sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )       
292         sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
293         IF( ln_zps ) THEN               ! z-coord. partial steps
294            DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
295               DO ji = 1, jpi
296                  ik = mbkt(ji,jj)
297                  IF( ik > 1 ) THEN
298                     zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
299                     sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
300                  ENDIF
301               END DO
302            END DO
303         ENDIF
304      ENDIF
305      !
306      CALL wrk_dealloc( jpi, jpj, jpk, 2, zsaldta )
307      !
308      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init')
309      !
310   END SUBROUTINE dia_ar5_init
311
312#else
313   !!----------------------------------------------------------------------
314   !!   Default option :                                         NO diaar5
315   !!----------------------------------------------------------------------
316   LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .FALSE.   ! coupled flag
317CONTAINS
318   SUBROUTINE dia_ar5_init    ! Dummy routine
319   END SUBROUTINE dia_ar5_init
320   SUBROUTINE dia_ar5( kt )   ! Empty routine
321      INTEGER ::   kt
322      WRITE(*,*) 'dia_ar5: You should not have seen this print! error?', kt
323   END SUBROUTINE dia_ar5
324#endif
325
326   !!======================================================================
327END MODULE diaar5
Note: See TracBrowser for help on using the repository browser.