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 NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA/diaar5.F90 @ 12344

Last change on this file since 12344 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

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