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.
diawri.F90 in NEMO/branches/UKMO/NEMO_4.0_surge/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_surge/src/OCE/DIA/diawri.F90 @ 11180

Last change on this file since 11180 was 11180, checked in by clne, 5 years ago

Initial commit of code for 2d (surge) work in NEMO4.
This is aiming to replicate the 3.6 version in branches/UKMO/dev_r5518_Surge_Modelling

File size: 49.7 KB
RevLine 
[3]1MODULE diawri
2   !!======================================================================
3   !!                     ***  MODULE  diawri  ***
4   !! Ocean diagnostics :  write ocean output files
5   !!=====================================================================
[2528]6   !! History :  OPA  ! 1991-03  (M.-A. Foujols)  Original code
7   !!            4.0  ! 1991-11  (G. Madec)
8   !!                 ! 1992-06  (M. Imbard)  correction restart file
9   !!                 ! 1992-07  (M. Imbard)  split into diawri and rstwri
10   !!                 ! 1993-03  (M. Imbard)  suppress writibm
11   !!                 ! 1998-01  (C. Levy)  NETCDF format using ioipsl INTERFACE
12   !!                 ! 1999-02  (E. Guilyardi)  name of netCDF files + variables
13   !!            8.2  ! 2000-06  (M. Imbard)  Original code (diabort.F)
14   !!   NEMO     1.0  ! 2002-06  (A.Bozec, E. Durand)  Original code (diainit.F)
15   !!             -   ! 2002-09  (G. Madec)  F90: Free form and module
16   !!             -   ! 2002-12  (G. Madec)  merge of diabort and diainit, F90
17   !!                 ! 2005-11  (V. Garnier) Surface pressure gradient organization
18   !!            3.2  ! 2008-11  (B. Lemaire) creation from old diawri
[5836]19   !!            3.7  ! 2014-01  (G. Madec) remove eddy induced velocity from no-IOM output
20   !!                 !                     change name of output variables in dia_wri_state
[2528]21   !!----------------------------------------------------------------------
[3]22
23   !!----------------------------------------------------------------------
[2528]24   !!   dia_wri       : create the standart output files
25   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
26   !!----------------------------------------------------------------------
[9019]27   USE oce            ! ocean dynamics and tracers
28   USE dom_oce        ! ocean space and time domain
29   USE phycst         ! physical constants
30   USE dianam         ! build name of file (routine)
31   USE diahth         ! thermocline diagnostics
32   USE dynadv   , ONLY: ln_dynadv_vec
33   USE icb_oce        ! Icebergs
34   USE icbdia         ! Iceberg budgets
35   USE ldftra         ! lateral physics: eddy diffusivity coef.
36   USE ldfdyn         ! lateral physics: eddy viscosity   coef.
37   USE sbc_oce        ! Surface boundary condition: ocean fields
38   USE sbc_ice        ! Surface boundary condition: ice fields
39   USE sbcssr         ! restoring term toward SST/SSS climatology
40   USE sbcwave        ! wave parameters
41   USE wet_dry        ! wetting and drying
42   USE zdf_oce        ! ocean vertical physics
43   USE zdfdrg         ! ocean vertical physics: top/bottom friction
44   USE zdfmxl         ! mixed layer
[6140]45   !
[9019]46   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
47   USE in_out_manager ! I/O manager
48   USE diatmb         ! Top,middle,bottom output
49   USE dia25h         ! 25h Mean output
50   USE iom            !
51   USE ioipsl         !
[5463]52
[9570]53#if defined key_si3
[10425]54   USE ice 
[9019]55   USE icewri 
[1482]56#endif
[2715]57   USE lib_mpp         ! MPP library
[3294]58   USE timing          ! preformance summary
[6140]59   USE diurnal_bulk    ! diurnal warm layer
60   USE cool_skin       ! Cool skin
[2528]61
[3]62   IMPLICIT NONE
63   PRIVATE
64
[2528]65   PUBLIC   dia_wri                 ! routines called by step.F90
66   PUBLIC   dia_wri_state
[2715]67   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
[3]68
[2528]69   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
[3609]70   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
[2528]71   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
72   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
73   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
74   INTEGER ::   ndex(1)                              ! ???
[2715]75   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
76   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
[3609]77   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
[3]78
79   !! * Substitutions
[1756]80#  include "vectopt_loop_substitute.h90"
[3]81   !!----------------------------------------------------------------------
[9598]82   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]83   !! $Id$
[10068]84   !! Software governed by the CeCILL license (see ./LICENSE)
[3]85   !!----------------------------------------------------------------------
86CONTAINS
87
[6140]88#if defined key_iomput
[3]89   !!----------------------------------------------------------------------
[2528]90   !!   'key_iomput'                                        use IOM library
91   !!----------------------------------------------------------------------
[9652]92   INTEGER FUNCTION dia_wri_alloc()
93      !
94      dia_wri_alloc = 0
95      !
96   END FUNCTION dia_wri_alloc
[2715]97
[9652]98   
[1561]99   SUBROUTINE dia_wri( kt )
[1482]100      !!---------------------------------------------------------------------
101      !!                  ***  ROUTINE dia_wri  ***
102      !!                   
103      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
104      !!      NETCDF format is used by default
105      !!
106      !! ** Method  :  use iom_put
107      !!----------------------------------------------------------------------
108      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
[1756]109      !!
[9019]110      INTEGER ::   ji, jj, jk       ! dummy loop indices
111      INTEGER ::   ikbot            ! local integer
112      REAL(wp)::   zztmp , zztmpx   ! local scalar
113      REAL(wp)::   zztmp2, zztmpy   !   -      -
114      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
115      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
[1482]116      !!----------------------------------------------------------------------
117      !
[9124]118      IF( ln_timing )   CALL timing_start('dia_wri')
[3294]119      !
[1482]120      ! Output the initial state and forcings
121      IF( ninist == 1 ) THEN                       
[10425]122         CALL dia_wri_state( 'output.init' )
[1482]123         ninist = 0
124      ENDIF
[3]125
[6351]126      ! Output of initial vertical scale factor
127      CALL iom_put("e3t_0", e3t_0(:,:,:) )
[10114]128      CALL iom_put("e3u_0", e3u_0(:,:,:) )
129      CALL iom_put("e3v_0", e3v_0(:,:,:) )
[6351]130      !
[6387]131      CALL iom_put( "e3t" , e3t_n(:,:,:) )
132      CALL iom_put( "e3u" , e3u_n(:,:,:) )
133      CALL iom_put( "e3v" , e3v_n(:,:,:) )
134      CALL iom_put( "e3w" , e3w_n(:,:,:) )
[6351]135      IF( iom_use("e3tdef") )   &
[6387]136         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
[5461]137
[9023]138      IF( ll_wd ) THEN
139         CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
140      ELSE
141         CALL iom_put( "ssh" , sshn )              ! sea surface height
142      ENDIF
143
[7646]144      IF( iom_use("wetdep") )   &                  ! wet depth
[9023]145         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
[5107]146     
147      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
148      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
149      IF ( iom_use("sbt") ) THEN
[4990]150         DO jj = 1, jpj
151            DO ji = 1, jpi
[9019]152               ikbot = mbkt(ji,jj)
153               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
[4990]154            END DO
[5107]155         END DO
156         CALL iom_put( "sbt", z2d )                ! bottom temperature
157      ENDIF
158     
159      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
160      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
161      IF ( iom_use("sbs") ) THEN
[4990]162         DO jj = 1, jpj
163            DO ji = 1, jpi
[9019]164               ikbot = mbkt(ji,jj)
165               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
[4990]166            END DO
[5107]167         END DO
168         CALL iom_put( "sbs", z2d )                ! bottom salinity
169      ENDIF
[5463]170
171      IF ( iom_use("taubot") ) THEN                ! bottom stress
[9019]172         zztmp = rau0 * 0.25
[7753]173         z2d(:,:) = 0._wp
[5463]174         DO jj = 2, jpjm1
175            DO ji = fs_2, fs_jpim1   ! vector opt.
[9019]176               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   &
177                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   &
178                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   &
179                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2
180               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
[5463]181               !
[9019]182            END DO
183         END DO
[10425]184         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
[5463]185         CALL iom_put( "taubot", z2d )           
186      ENDIF
[11180]187 
188      IF( iom_use("uwnd") ) CALL iom_put( "uwnd" ,   uwnd  ) 
189      IF( iom_use("vwnd") ) CALL iom_put( "vwnd" ,   vwnd  ) 
190         
[9019]191      CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current
192      CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current
[5107]193      IF ( iom_use("sbu") ) THEN
[4990]194         DO jj = 1, jpj
195            DO ji = 1, jpi
[9019]196               ikbot = mbku(ji,jj)
197               z2d(ji,jj) = un(ji,jj,ikbot)
[4990]198            END DO
[5107]199         END DO
200         CALL iom_put( "sbu", z2d )                ! bottom i-current
201      ENDIF
202     
[9019]203      CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current
204      CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current
[5107]205      IF ( iom_use("sbv") ) THEN
[4990]206         DO jj = 1, jpj
207            DO ji = 1, jpi
[9019]208               ikbot = mbkv(ji,jj)
209               z2d(ji,jj) = vn(ji,jj,ikbot)
[4990]210            END DO
[5107]211         END DO
212         CALL iom_put( "sbv", z2d )                ! bottom j-current
[4990]213      ENDIF
[1482]214
[5461]215      CALL iom_put( "woce", wn )                   ! vertical velocity
216      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
217         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
[7753]218         z2d(:,:) = rau0 * e1e2t(:,:)
[5461]219         DO jk = 1, jpk
[7753]220            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
[5461]221         END DO
222         CALL iom_put( "w_masstr" , z3d ) 
223         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
224      ENDIF
225
[9019]226      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
227      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef.
228      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef.
[5107]229
[9019]230      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
231      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
[6351]232
[5107]233      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
[4990]234         DO jj = 2, jpjm1                                    ! sst gradient
235            DO ji = fs_2, fs_jpim1   ! vector opt.
[5836]236               zztmp  = tsn(ji,jj,1,jp_tem)
237               zztmpx = ( tsn(ji+1,jj,1,jp_tem) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) * r1_e1u(ji-1,jj)
238               zztmpy = ( tsn(ji,jj+1,1,jp_tem) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) * r1_e2v(ji,jj-1)
[4990]239               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
240                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
241            END DO
[1756]242         END DO
[10425]243         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
[9019]244         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient
[7753]245         z2d(:,:) = SQRT( z2d(:,:) )
[9019]246         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient
[4990]247      ENDIF
248         
[9019]249      ! heat and salt contents
[4990]250      IF( iom_use("heatc") ) THEN
[7753]251         z2d(:,:)  = 0._wp 
[4990]252         DO jk = 1, jpkm1
[5107]253            DO jj = 1, jpj
254               DO ji = 1, jpi
[6140]255                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
[4990]256               END DO
[4761]257            END DO
258         END DO
[9019]259         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
[4990]260      ENDIF
261
262      IF( iom_use("saltc") ) THEN
[7753]263         z2d(:,:)  = 0._wp 
[4990]264         DO jk = 1, jpkm1
[5107]265            DO jj = 1, jpj
266               DO ji = 1, jpi
[6140]267                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
[4990]268               END DO
269            END DO
270         END DO
[9019]271         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
[4990]272      ENDIF
[4840]273      !
[4990]274      IF ( iom_use("eken") ) THEN
[9399]275         z3d(:,:,jpk) = 0._wp 
[4990]276         DO jk = 1, jpkm1
277            DO jj = 2, jpjm1
278               DO ji = fs_2, fs_jpim1   ! vector opt.
[9019]279                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
280                  z3d(ji,jj,jk) = zztmp * (  un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   &
281                     &                     + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   &
282                     &                     + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   &
283                     &                     + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   )
284               END DO
285            END DO
286         END DO
[10425]287         CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
[9019]288         CALL iom_put( "eken", z3d )                 ! kinetic energy
[4990]289      ENDIF
[6351]290      !
291      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
292      !
[7646]293      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
[7753]294         z3d(:,:,jpk) = 0.e0
295         z2d(:,:) = 0.e0
[1756]296         DO jk = 1, jpkm1
[7753]297            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
298            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
[1756]299         END DO
[9019]300         CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction
301         CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum
[4990]302      ENDIF
303     
304      IF( iom_use("u_heattr") ) THEN
[9019]305         z2d(:,:) = 0._wp 
[4990]306         DO jk = 1, jpkm1
307            DO jj = 2, jpjm1
308               DO ji = fs_2, fs_jpim1   ! vector opt.
309                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
310               END DO
311            END DO
312         END DO
[10425]313         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
[9019]314         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction
[4990]315      ENDIF
[4761]316
[4990]317      IF( iom_use("u_salttr") ) THEN
[7753]318         z2d(:,:) = 0.e0 
[1756]319         DO jk = 1, jpkm1
320            DO jj = 2, jpjm1
321               DO ji = fs_2, fs_jpim1   ! vector opt.
[4990]322                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
[1756]323               END DO
324            END DO
325         END DO
[10425]326         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
[9019]327         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction
[4990]328      ENDIF
[4761]329
[4990]330     
331      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
[7753]332         z3d(:,:,jpk) = 0.e0
[1756]333         DO jk = 1, jpkm1
[7753]334            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
[1756]335         END DO
[9019]336         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
[4990]337      ENDIF
338     
339      IF( iom_use("v_heattr") ) THEN
[7753]340         z2d(:,:) = 0.e0 
[4990]341         DO jk = 1, jpkm1
342            DO jj = 2, jpjm1
343               DO ji = fs_2, fs_jpim1   ! vector opt.
344                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
345               END DO
346            END DO
347         END DO
[10425]348         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
[9019]349         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction
[4990]350      ENDIF
[4761]351
[4990]352      IF( iom_use("v_salttr") ) THEN
[9019]353         z2d(:,:) = 0._wp 
[1756]354         DO jk = 1, jpkm1
355            DO jj = 2, jpjm1
356               DO ji = fs_2, fs_jpim1   ! vector opt.
[4990]357                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
[1756]358               END DO
359            END DO
360         END DO
[10425]361         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
[9019]362         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction
[1756]363      ENDIF
[7646]364
365      IF( iom_use("tosmint") ) THEN
[9019]366         z2d(:,:) = 0._wp
[7646]367         DO jk = 1, jpkm1
368            DO jj = 2, jpjm1
369               DO ji = fs_2, fs_jpim1   ! vector opt.
[9019]370                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
[7646]371               END DO
372            END DO
373         END DO
[10425]374         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
[9019]375         CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature
[7646]376      ENDIF
377      IF( iom_use("somint") ) THEN
[7753]378         z2d(:,:)=0._wp
[7646]379         DO jk = 1, jpkm1
380            DO jj = 2, jpjm1
381               DO ji = fs_2, fs_jpim1   ! vector opt.
[9019]382                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
[7646]383               END DO
384            END DO
385         END DO
[10425]386         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
[9019]387         CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity
[7646]388      ENDIF
389
[9019]390      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2)
[2528]391      !
[6140]392
[9019]393      IF (ln_diatmb)   CALL dia_tmb                   ! tmb values
394         
395      IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging
[6140]396
[9124]397      IF( ln_timing )   CALL timing_stop('dia_wri')
[3294]398      !
[1482]399   END SUBROUTINE dia_wri
400
401#else
[2528]402   !!----------------------------------------------------------------------
403   !!   Default option                                  use IOIPSL  library
404   !!----------------------------------------------------------------------
405
[9652]406   INTEGER FUNCTION dia_wri_alloc()
407      !!----------------------------------------------------------------------
408      INTEGER, DIMENSION(2) :: ierr
409      !!----------------------------------------------------------------------
410      ierr = 0
411      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
412         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
413         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
414         !
415      dia_wri_alloc = MAXVAL(ierr)
[10425]416      CALL mpp_sum( 'diawri', dia_wri_alloc )
[9652]417      !
418   END FUNCTION dia_wri_alloc
419
420   
[1561]421   SUBROUTINE dia_wri( kt )
[3]422      !!---------------------------------------------------------------------
423      !!                  ***  ROUTINE dia_wri  ***
424      !!                   
425      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
426      !!      NETCDF format is used by default
427      !!
428      !! ** Method  :   At the beginning of the first time step (nit000),
429      !!      define all the NETCDF files and fields
430      !!      At each time step call histdef to compute the mean if ncessary
431      !!      Each nwrite time step, output the instantaneous or mean fields
432      !!----------------------------------------------------------------------
[5836]433      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
434      !
[2528]435      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
436      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
437      INTEGER  ::   inum = 11                                ! temporary logical unit
[3294]438      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
439      INTEGER  ::   ierr                                     ! error code return from allocation
[2528]440      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
[3609]441      INTEGER  ::   jn, ierror                               ! local integers
[6140]442      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
[5836]443      !
[9019]444      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
445      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
[3]446      !!----------------------------------------------------------------------
[3294]447      !
[9124]448      IF( ln_timing )   CALL timing_start('dia_wri')
[1482]449      !
[9019]450      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
[10425]451         CALL dia_wri_state( 'output.init' )
[1482]452         ninist = 0
453      ENDIF
454      !
[3]455      ! 0. Initialisation
456      ! -----------------
[632]457
[9019]458      ll_print = .FALSE.                  ! local variable for debugging
[3]459      ll_print = ll_print .AND. lwp
460
461      ! Define frequency of output and means
[5566]462      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
[3]463#if defined key_diainstant
[6140]464      zsto = nwrite * rdt
[1312]465      clop = "inst("//TRIM(clop)//")"
[3]466#else
[6140]467      zsto=rdt
[1312]468      clop = "ave("//TRIM(clop)//")"
[3]469#endif
[6140]470      zout = nwrite * rdt
471      zmax = ( nitend - nit000 + 1 ) * rdt
[3]472
473      ! Define indices of the horizontal output zoom and vertical limit storage
474      iimi = 1      ;      iima = jpi
475      ijmi = 1      ;      ijma = jpj
476      ipk = jpk
477
478      ! define time axis
[1334]479      it = kt
480      itmod = kt - nit000 + 1
[3]481
482
483      ! 1. Define NETCDF files and fields at beginning of first time step
484      ! -----------------------------------------------------------------
485
486      IF( kt == nit000 ) THEN
487
488         ! Define the NETCDF files (one per grid)
[632]489
[3]490         ! Compute julian date from starting date of the run
[1309]491         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
492         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
[3]493         IF(lwp)WRITE(numout,*)
494         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
495            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
496         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
497                                 ' limit storage in depth = ', ipk
498
499         ! WRITE root name in date.file for use by postpro
[1581]500         IF(lwp) THEN
[895]501            CALL dia_nam( clhstnam, nwrite,' ' )
[1581]502            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]503            WRITE(inum,*) clhstnam
504            CLOSE(inum)
505         ENDIF
[632]506
[3]507         ! Define the T grid FILE ( nid_T )
[632]508
[3]509         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
510         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
511         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
512            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
[6140]513            &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]514         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
[4292]515            &           "m", ipk, gdept_1d, nz_T, "down" )
[3]516         !                                                            ! Index of ocean points
517         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
518         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
[3609]519         !
520         IF( ln_icebergs ) THEN
521            !
522            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
523            !! that routine is called from nemogcm, so do it here immediately before its needed
524            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
[10425]525            CALL mpp_sum( 'diawri', ierror )
[3609]526            IF( ierror /= 0 ) THEN
527               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
528               RETURN
529            ENDIF
530            !
531            !! iceberg vertical coordinate is class number
532            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
533               &           "number", nclasses, class_num, nb_T )
534            !
535            !! each class just needs the surface index pattern
536            ndim_bT = 3
537            DO jn = 1,nclasses
538               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
539            ENDDO
540            !
541         ENDIF
[3]542
543         ! Define the U grid FILE ( nid_U )
544
545         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
546         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
547         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
548            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
[6140]549            &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]550         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
[4292]551            &           "m", ipk, gdept_1d, nz_U, "down" )
[3]552         !                                                            ! Index of ocean points
553         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
554         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
555
556         ! Define the V grid FILE ( nid_V )
557
558         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
559         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
560         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
561            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
[6140]562            &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]563         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
[4292]564            &          "m", ipk, gdept_1d, nz_V, "down" )
[3]565         !                                                            ! Index of ocean points
566         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
567         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
568
569         ! Define the W grid FILE ( nid_W )
570
571         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
572         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
573         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
574            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
[6140]575            &          nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]576         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
[4292]577            &          "m", ipk, gdepw_1d, nz_W, "down" )
[3]578
[632]579
[3]580         ! Declare all the output fields as NETCDF variables
581
582         !                                                                                      !!! nid_T : 3D
583         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
584            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
585         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
586            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[6140]587         IF(  .NOT.ln_linssh  ) THEN
[4292]588            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
589            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
590            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
591            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
592            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
593            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
594         ENDIF
[3]595         !                                                                                      !!! nid_T : 2D
596         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
597            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
598         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
599            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
[359]600         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
[3]601            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
[2528]602         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
[3]603            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
[4570]604         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
605            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
[3625]606         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
[3]607            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
[6140]608         IF(  ln_linssh  ) THEN
[4292]609            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
[3625]610            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
[4292]611            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
612            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
[3625]613            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
[4292]614            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
615         ENDIF
[888]616         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
[3]617            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
618         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
619            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
[1585]620         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
621            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
[3]622         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
623            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
[1037]624         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
[3]625            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
[1649]626         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
627            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
[3609]628!
629         IF( ln_icebergs ) THEN
630            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
631               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
632            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
633               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
634            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
635               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
636            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
637               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
638            IF( ln_bergdia ) THEN
639               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
640                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
641               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
642                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
643               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
644                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
645               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
646                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
647               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
648                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
649               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
650                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
651               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
652                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
653               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
654                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
655               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
656                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
657               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
658                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
659            ENDIF
660         ENDIF
661
[5407]662         IF( .NOT. ln_cpl ) THEN
[4990]663            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
664               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
665            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
666               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
667            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
668               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
669         ENDIF
[3]670
[5407]671         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
[4990]672            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
673               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
674            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
675               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
676            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
677               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
678         ENDIF
679         
[3]680         clmx ="l_max(only(x))"    ! max index on a period
[5836]681!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
682!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
[3]683#if defined key_diahth
684         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
685            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
686         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
687            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
688         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
689            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
[7646]690         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "J/m2"   ,   & ! htc3
[3]691            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
692#endif
693
[2528]694         CALL histend( nid_T, snc4chunks=snc4set )
[3]695
696         !                                                                                      !!! nid_U : 3D
697         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
698            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[7646]699         IF( ln_wave .AND. ln_sdw) THEN
700            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd
701               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
702         ENDIF
[3]703         !                                                                                      !!! nid_U : 2D
[888]704         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
[3]705            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
706
[2528]707         CALL histend( nid_U, snc4chunks=snc4set )
[3]708
709         !                                                                                      !!! nid_V : 3D
710         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
711            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[7646]712         IF( ln_wave .AND. ln_sdw) THEN
713            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd
714               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
715         ENDIF
[3]716         !                                                                                      !!! nid_V : 2D
[888]717         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
[3]718            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
719
[2528]720         CALL histend( nid_V, snc4chunks=snc4set )
[3]721
722         !                                                                                      !!! nid_W : 3D
723         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
724            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
725         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
726            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[9019]727         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avm
[255]728            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
729
[9019]730         IF( ln_zdfddm ) THEN
[3]731            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
732               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
733         ENDIF
[7646]734         
735         IF( ln_wave .AND. ln_sdw) THEN
736            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd
737               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
738         ENDIF
[3]739         !                                                                                      !!! nid_W : 2D
[2528]740         CALL histend( nid_W, snc4chunks=snc4set )
[3]741
742         IF(lwp) WRITE(numout,*)
743         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
744         IF(ll_print) CALL FLUSH(numout )
745
746      ENDIF
747
748      ! 2. Start writing data
749      ! ---------------------
750
[4292]751      ! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]752      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
753      ! donne le nombre d'elements, et ndex la liste des indices a sortir
754
[1334]755      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
[3]756         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
757         WRITE(numout,*) '~~~~~~ '
758      ENDIF
759
[6140]760      IF( .NOT.ln_linssh ) THEN
761         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
762         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
763         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
764         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
[4292]765      ELSE
766         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
767         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
768         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
769         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
770      ENDIF
[6140]771      IF( .NOT.ln_linssh ) THEN
[7753]772         zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
[6140]773         CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
774         CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
[4292]775         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
776      ENDIF
[3]777      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
[2528]778      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
[4570]779      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
[3625]780      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
781                                                                                  ! (includes virtual salt flux beneath ice
782                                                                                  ! in linear free surface case)
[6140]783      IF( ln_linssh ) THEN
[7753]784         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
[4292]785         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
[7753]786         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
[4292]787         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
788      ENDIF
[888]789      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
[3]790      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
[1585]791      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
[3]792      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
[1037]793      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
[1649]794      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
[3609]795!
796      IF( ln_icebergs ) THEN
797         !
798         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
799         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
800         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
801         !
802         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
803         !
804         IF( ln_bergdia ) THEN
805            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
806            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
807            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
808            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
809            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
810            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
811            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
812            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
813            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
814            !
815            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
816         ENDIF
817      ENDIF
818
[5407]819      IF( .NOT. ln_cpl ) THEN
[4990]820         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
821         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
[7753]822         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[4990]823         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
824      ENDIF
[5407]825      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
[4990]826         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
827         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
[7753]828         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[4990]829         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
830      ENDIF
831!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
832!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
[3]833
834#if defined key_diahth
835      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
836      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
837      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
838      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
839#endif
[888]840
[3]841      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
[888]842      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
[3]843
844      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
[888]845      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
[3]846
847      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
848      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
[9019]849      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
850      IF( ln_zdfddm ) THEN
851         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef.
[3]852      ENDIF
853
[7646]854      IF( ln_wave .AND. ln_sdw ) THEN
[9019]855         CALL histwrite( nid_U, "sdzocrtx", it, usd         , ndim_U , ndex_U )    ! i-StokesDrift-current
856         CALL histwrite( nid_V, "sdmecrty", it, vsd         , ndim_V , ndex_V )    ! j-StokesDrift-current
857         CALL histwrite( nid_W, "sdvecrtz", it, wsd         , ndim_T , ndex_T )    ! StokesDrift vert. current
[7646]858      ENDIF
859
[1318]860      ! 3. Close all files
861      ! ---------------------------------------
[1561]862      IF( kt == nitend ) THEN
[3]863         CALL histclo( nid_T )
864         CALL histclo( nid_U )
865         CALL histclo( nid_V )
866         CALL histclo( nid_W )
867      ENDIF
[2528]868      !
[9124]869      IF( ln_timing )   CALL timing_stop('dia_wri')
[3294]870      !
[3]871   END SUBROUTINE dia_wri
[1567]872#endif
873
[10425]874   SUBROUTINE dia_wri_state( cdfile_name )
[3]875      !!---------------------------------------------------------------------
876      !!                 ***  ROUTINE dia_wri_state  ***
877      !!       
878      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
879      !!      the instantaneous ocean state and forcing fields.
880      !!        Used to find errors in the initial state or save the last
881      !!      ocean state in case of abnormal end of a simulation
882      !!
883      !! ** Method  :   NetCDF files using ioipsl
884      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
885      !!      File 'output.abort.nc' is created in case of abnormal job end
886      !!----------------------------------------------------------------------
[1334]887      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
[10425]888      !!
889      INTEGER :: inum
[3]890      !!----------------------------------------------------------------------
[3294]891      !
[648]892      IF(lwp) WRITE(numout,*)
893      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
894      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
[10425]895      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
[648]896
[9570]897#if defined key_si3
[10425]898     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
[1482]899#else
[10425]900     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
[1482]901#endif
[3]902
[10425]903      CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature
904      CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity
905      CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height
906      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
907      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
908      CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity
[7646]909      IF( ALLOCATED(ahtu) ) THEN
[10425]910         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
911         CALL iom_rstput( 0, 0, inum,  'ahtv', ahtv              )    ! aht at v-point
[7646]912      ENDIF
913      IF( ALLOCATED(ahmt) ) THEN
[10425]914         CALL iom_rstput( 0, 0, inum,  'ahmt', ahmt              )    ! ahmt at u-point
915         CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
[7646]916      ENDIF
[10425]917      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
918      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
919      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
920      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
921      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
922      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
[6140]923      IF(  .NOT.ln_linssh  ) THEN             
[10425]924         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n        )    !  T-cell depth
925         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n          )    !  T-cell thickness 
926      END IF
[7646]927      IF( ln_wave .AND. ln_sdw ) THEN
[10425]928         CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
929         CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd            )    ! now StokesDrift j-velocity
930         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
[7646]931      ENDIF
[10425]932 
933#if defined key_si3
934      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
935         CALL ice_wri_state( inum )
[1561]936      ENDIF
937#endif
[10425]938      !
939      CALL iom_close( inum )
[3294]940      !
[3]941   END SUBROUTINE dia_wri_state
[6140]942
[3]943   !!======================================================================
944END MODULE diawri
Note: See TracBrowser for help on using the repository browser.