source: NEMO/trunk/src/OCE/DIA/diaar5.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 10 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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