source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 15.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   !!   dia_ar5       : AR5 diagnostics
10   !!   dia_ar5_init  : initialisation of AR5 diagnostics
11   !!----------------------------------------------------------------------
12   USE oce            ! ocean dynamics and active tracers
13   USE dom_oce        ! ocean space and time domain
14   USE eosbn2         ! equation of state                (eos_bn2 routine)
15   USE lib_mpp        ! distribued memory computing library
16   USE iom            ! I/O manager library
17   USE timing         ! preformance summary
18   USE fldread        ! type FLD_N
19   USE phycst         ! physical constant
20   USE in_out_manager  ! I/O manager
21   USE zdfddm
22   USE zdf_oce
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   dia_ar5        ! routine called in step.F90 module
28   PUBLIC   dia_ar5_alloc  ! routine called in nemogcm.F90 module
29   PUBLIC   dia_ar5_hst    ! heat/salt transport
30
31   REAL(wp)                         ::   vol0         ! ocean volume (interior domain)
32   REAL(wp)                         ::   area_tot     ! total ocean surface (interior domain)
33   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   area         ! cell surface (interior domain)
34   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   thick0       ! ocean thickness (interior domain)
35   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity
36
37   LOGICAL  :: l_ar5
38     
39   !! * Substitutions
40#  include "zdfddm_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   FUNCTION dia_ar5_alloc()
50      !!----------------------------------------------------------------------
51      !!                    ***  ROUTINE dia_ar5_alloc  ***
52      !!----------------------------------------------------------------------
53      INTEGER :: dia_ar5_alloc
54      !!----------------------------------------------------------------------
55      !
56      ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc )
57      !
58      IF( lk_mpp             )   CALL mpp_sum ( dia_ar5_alloc )
59      IF( dia_ar5_alloc /= 0 )   CALL ctl_warn('dia_ar5_alloc: failed to allocate arrays')
60      !
61   END FUNCTION dia_ar5_alloc
62
63
64   SUBROUTINE dia_ar5( kt )
65      !!----------------------------------------------------------------------
66      !!                    ***  ROUTINE dia_ar5  ***
67      !!
68      !! ** Purpose :   compute and output some AR5 diagnostics
69      !!----------------------------------------------------------------------
70      !
71      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
72      !
73      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
74      REAL(wp) ::   zvolssh, zvol, zssh_steric, zztmp, zarho, ztemp, zsal, zmass
75      REAL(wp) ::   zaw, zbw, zrw
76      !
77      REAL(wp), DIMENSION(jpi,jpj)     :: zarea_ssh , zbotpres       ! 2D workspace
78      REAL(wp), DIMENSION(jpi,jpj)     :: zpe                         ! 2D workspace
79      REAL(wp), DIMENSION(jpi,jpj,jpk)   :: zrhd , zrhop               ! 3D workspace
80      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: ztsn                       ! 4D workspace
81      !!--------------------------------------------------------------------
82      IF( nn_timing == 1 )   CALL timing_start('dia_ar5')
83 
84      IF( kt == nit000 )     CALL dia_ar5_init
85
86      IF( l_ar5 ) THEN
87         zarea_ssh(:,:) = area(:,:) * sshn(:,:)
88      ENDIF
89      !
90      IF( iom_use( 'voltot' ) .OR. iom_use( 'sshtot' )  .OR. iom_use( 'sshdyn' )  ) THEN   
91         !                                         ! total volume of liquid seawater
92         zvolssh = SUM( zarea_ssh(:,:) ) 
93         IF( lk_mpp )   CALL mpp_sum( zvolssh )
94         zvol = vol0 + zvolssh
95     
96         CALL iom_put( 'voltot', zvol               )
97         CALL iom_put( 'sshtot', zvolssh / area_tot )
98         CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) )
99         !
100      ENDIF
101
102      IF( iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) THEN   
103         !                     
104         ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh
105         ztsn(:,:,:,jp_sal) = sn0(:,:,:)
106         CALL eos( ztsn, zrhd, gdept_n(:,:,:) )                       ! now in situ density using initial salinity
107         !
108         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
109         DO jk = 1, jpkm1
110            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
111         END DO
112         IF( ln_linssh ) THEN
113            IF( ln_isfcav ) THEN
114               DO ji = 1, jpi
115                  DO jj = 1, jpj
116                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
117                  END DO
118               END DO
119            ELSE
120               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
121            END IF
122!!gm
123!!gm   riceload should be added in both ln_linssh=T or F, no?
124!!gm
125         END IF
126         !                                         
127         zarho = SUM( area(:,:) * zbotpres(:,:) ) 
128         IF( lk_mpp )   CALL mpp_sum( zarho )
129         zssh_steric = - zarho / area_tot
130         CALL iom_put( 'sshthster', zssh_steric )
131     
132         !                                         ! steric sea surface height
133         CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density
134         zrhop(:,:,jpk) = 0._wp
135         CALL iom_put( 'rhop', zrhop )
136         !
137         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
138         DO jk = 1, jpkm1
139            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
140         END DO
141         IF( ln_linssh ) THEN
142            IF ( ln_isfcav ) THEN
143               DO ji = 1,jpi
144                  DO jj = 1,jpj
145                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
146                  END DO
147               END DO
148            ELSE
149               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
150            END IF
151         END IF
152         !   
153         zarho = SUM( area(:,:) * zbotpres(:,:) ) 
154         IF( lk_mpp )   CALL mpp_sum( zarho )
155         zssh_steric = - zarho / area_tot
156         CALL iom_put( 'sshsteric', zssh_steric )
157     
158         !                                         ! ocean bottom pressure
159         zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
160         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
161         CALL iom_put( 'botpres', zbotpres )
162         !
163      ENDIF
164
165      IF( iom_use( 'masstot' ) .OR. iom_use( 'temptot' )  .OR. iom_use( 'saltot' )  ) THEN   
166         !                                         ! Mean density anomalie, temperature and salinity
167         ztemp = 0._wp
168         zsal  = 0._wp
169         DO jk = 1, jpkm1
170            DO jj = 1, jpj
171               DO ji = 1, jpi
172                  zztmp = area(ji,jj) * e3t_n(ji,jj,jk)
173                  ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem)
174                  zsal  = zsal  + zztmp * tsn(ji,jj,jk,jp_sal)
175               END DO
176            END DO
177         END DO
178         IF( ln_linssh ) THEN
179            IF( ln_isfcav ) THEN
180               DO ji = 1, jpi
181                  DO jj = 1, jpj
182                     ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem) 
183                     zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal) 
184                  END DO
185               END DO
186            ELSE
187               ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
188               zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
189            END IF
190         ENDIF
191         IF( lk_mpp ) THEN 
192            CALL mpp_sum( ztemp )
193            CALL mpp_sum( zsal  )
194         END IF
195         !
196         zmass = rau0 * ( zarho + zvol )                 ! total mass of liquid seawater
197         ztemp = ztemp / zvol                            ! potential temperature in liquid seawater
198         zsal  = zsal  / zvol                            ! Salinity of liquid seawater
199         !
200         CALL iom_put( 'masstot', zmass )
201         CALL iom_put( 'temptot', ztemp )
202         CALL iom_put( 'saltot' , zsal  )
203         !
204      ENDIF
205
206      IF( iom_use( 'tnpeo' )) THEN   
207      ! Work done against stratification by vertical mixing
208      ! Exclude points where rn2 is negative as convection kicks in here and
209      ! work is not being done against stratification
210          zpe(:,:) = 0._wp
211          IF( lk_zdfddm ) THEN
212             DO ji=1,jpi
213                DO jj=1,jpj
214                   DO jk=1,jpk
215                      zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
216                         &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )
217                      !
218                      zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
219                      zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
220                      !
221                      zpe(ji, jj) = zpe(ji, jj) - MIN(0._wp, rn2(ji,jj,jk)) * &
222                           &       grav * (avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  &
223                           &       - fsavs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )
224
225                   ENDDO
226                ENDDO
227             ENDDO
228          ELSE
229             DO ji = 1, jpi
230                DO jj = 1, jpj
231                   DO jk = 1, jpk
232                       zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
233                   ENDDO
234                ENDDO
235             ENDDO
236          ENDIF
237          CALL lbc_lnk( zpe, 'T', 1._wp)         
238          CALL iom_put( 'tnpeo', zpe )
239      ENDIF
240      !
241      IF( l_ar5 ) THEN
242      ENDIF
243      !
244      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5')
245      !
246   END SUBROUTINE dia_ar5
247
248   SUBROUTINE dia_ar5_hst( ktra, cptr, pua, pva ) 
249      !!----------------------------------------------------------------------
250      !!                    ***  ROUTINE dia_ar5_htr ***
251      !!----------------------------------------------------------------------
252      !! Wrapper for heat transport calculations
253      !! Called from all advection and/or diffusion routines
254      !!----------------------------------------------------------------------
255      INTEGER                         , INTENT(in )  :: ktra  ! tracer index
256      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'
257      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pua   ! 3D input array of advection/diffusion
258      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion
259      !
260      INTEGER    ::  ji, jj, jk
261      REAL(wp), DIMENSION(jpi,jpj)  :: z2d
262
263   
264
265      z2d(:,:) = pua(:,:,1) 
266      DO jk = 1, jpkm1
267         DO jj = 2, jpjm1
268            DO ji = fs_2, fs_jpim1   ! vector opt.
269               z2d(ji,jj) = z2d(ji,jj) + pua(ji,jj,jk) 
270            END DO
271         END DO
272       END DO
273       CALL lbc_lnk( z2d, 'U', -1. )
274       IF( cptr == 'adv' ) THEN
275          IF( ktra == jp_tem ) CALL iom_put( "uadv_heattr" , rau0_rcp * z2d )  ! advective heat transport in i-direction
276          IF( ktra == jp_sal ) CALL iom_put( "uadv_salttr" , rau0     * z2d )  ! advective salt transport in i-direction
277       ENDIF
278       IF( cptr == 'ldf' ) THEN
279          IF( ktra == jp_tem ) CALL iom_put( "udiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in i-direction
280          IF( ktra == jp_sal ) CALL iom_put( "udiff_salttr" , rau0     * z2d ) ! diffusive salt transport in i-direction
281       ENDIF
282       !
283       z2d(:,:) = pva(:,:,1) 
284       DO jk = 1, jpkm1
285          DO jj = 2, jpjm1
286             DO ji = fs_2, fs_jpim1   ! vector opt.
287                z2d(ji,jj) = z2d(ji,jj) + pva(ji,jj,jk) 
288             END DO
289          END DO
290       END DO
291       CALL lbc_lnk( z2d, 'V', -1. )
292       IF( cptr == 'adv' ) THEN
293          IF( ktra == jp_tem ) CALL iom_put( "vadv_heattr" , rau0_rcp * z2d )  ! advective heat transport in j-direction
294          IF( ktra == jp_sal ) CALL iom_put( "vadv_salttr" , rau0     * z2d )  ! advective salt transport in j-direction
295       ENDIF
296       IF( cptr == 'ldf' ) THEN
297          IF( ktra == jp_tem ) CALL iom_put( "vdiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in j-direction
298          IF( ktra == jp_sal ) CALL iom_put( "vdiff_salttr" , rau0     * z2d ) ! diffusive salt transport in j-direction
299       ENDIF
300         
301
302   END SUBROUTINE dia_ar5_hst
303
304
305   SUBROUTINE dia_ar5_init
306      !!----------------------------------------------------------------------
307      !!                  ***  ROUTINE dia_ar5_init  ***
308      !!                   
309      !! ** Purpose :   initialization for AR5 diagnostic computation
310      !!----------------------------------------------------------------------
311      INTEGER  ::   inum
312      INTEGER  ::   ik
313      INTEGER  ::   ji, jj, jk  ! dummy loop indices
314      REAL(wp) ::   zztmp 
315      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::   zsaldta   ! Jan/Dec levitus salinity
316      !
317      !!----------------------------------------------------------------------
318      !
319      IF( nn_timing == 1 )   CALL timing_start('dia_ar5_init')
320      !
321      l_ar5 = .FALSE.
322      IF(   iom_use( 'voltot'  ) .OR. iom_use( 'sshtot'    )  .OR. iom_use( 'sshdyn' )  .OR.  & 
323         &  iom_use( 'masstot' ) .OR. iom_use( 'temptot'   )  .OR. iom_use( 'saltot' ) .OR.  &   
324         &  iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) L_ar5 = .TRUE.
325 
326      IF( l_ar5 ) THEN
327         !
328         !                                      ! allocate dia_ar5 arrays
329         IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
330
331         area(:,:) = e1e2t(:,:) * tmask_i(:,:)
332
333         area_tot = SUM( area(:,:) )   ;   IF( lk_mpp )   CALL mpp_sum( area_tot )
334
335         vol0        = 0._wp
336         thick0(:,:) = 0._wp
337         DO jk = 1, jpkm1
338            vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) )
339            thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk)
340         END DO
341         IF( lk_mpp )   CALL mpp_sum( vol0 )
342
343         CALL iom_open ( 'sali_ref_clim_monthly', inum )
344         CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  )
345         CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 )
346         CALL iom_close( inum )
347
348         sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )       
349         sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
350         IF( ln_zps ) THEN               ! z-coord. partial steps
351            DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
352               DO ji = 1, jpi
353                  ik = mbkt(ji,jj)
354                  IF( ik > 1 ) THEN
355                     zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
356                     sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
357                  ENDIF
358               END DO
359            END DO
360         ENDIF
361         !
362         !
363      ENDIF
364      !
365      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init')
366      !
367   END SUBROUTINE dia_ar5_init
368
369   !!======================================================================
370END MODULE diaar5
Note: See TracBrowser for help on using the repository browser.