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 branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 12963

Last change on this file since 12963 was 12963, checked in by clne, 4 years ago

Mask land in surge output fields with fill value rather than 0
  • Property svn:keywords set to Id
File size: 56.4 KB
Line 
1MODULE diawri
2   !!======================================================================
3   !!                     ***  MODULE  diawri  ***
4   !! Ocean diagnostics :  write ocean output files
5   !!=====================================================================
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
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   dia_wri       : create the standart output files
23   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE dynadv, ONLY: ln_dynadv_vec
28   USE zdf_oce         ! ocean vertical physics
29   USE ldftra_oce      ! ocean active tracers: lateral physics
30   USE ldfdyn_oce      ! ocean dynamics: lateral physics
31   USE traldf_iso_grif, ONLY : psix_eiv, psiy_eiv
32   USE sol_oce         ! solver variables
33   USE sbc_oce         ! Surface boundary condition: ocean fields
34   USE sbc_ice         ! Surface boundary condition: ice fields
35   USE icb_oce         ! Icebergs
36   USE icbdia          ! Iceberg budgets
37   USE sbcssr          ! restoring term toward SST/SSS climatology
38   USE phycst          ! physical constants
39   USE zdfmxl          ! mixed layer
40   USE dianam          ! build name of file (routine)
41   USE zdfddm          ! vertical  physics: double diffusion
42   USE diahth          ! thermocline diagnostics
43   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE in_out_manager  ! I/O manager
45   USE diadimg         ! dimg direct access file format output
46   USE iom
47   USE ioipsl
48   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities     
49
50#if defined key_lim2
51   USE limwri_2 
52#elif defined key_lim3
53   USE limwri 
54#endif
55   USE lib_mpp         ! MPP library
56   USE timing          ! preformance summary
57   USE wrk_nemo        ! working array
58
59   IMPLICIT NONE
60   PRIVATE
61
62   PUBLIC   dia_wri                 ! routines called by step.F90
63   PUBLIC   dia_wri_state
64   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
65
66   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
67   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
68   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
69   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
70   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
71   INTEGER ::   ndex(1)                              ! ???
72   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
73   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
74   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
75
76   !! * Substitutions
77#  include "zdfddm_substitute.h90"
78#  include "domzgr_substitute.h90"
79#  include "vectopt_loop_substitute.h90"
80   !!----------------------------------------------------------------------
81   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
82   !! $Id$
83   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
84   !!----------------------------------------------------------------------
85CONTAINS
86
87   INTEGER FUNCTION dia_wri_alloc()
88      !!----------------------------------------------------------------------
89      INTEGER, DIMENSION(2) :: ierr
90      !!----------------------------------------------------------------------
91      ierr = 0
92      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
93         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
94         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
95         !
96      dia_wri_alloc = MAXVAL(ierr)
97      IF( lk_mpp )   CALL mpp_sum( dia_wri_alloc )
98      !
99  END FUNCTION dia_wri_alloc
100
101#if defined key_dimgout
102   !!----------------------------------------------------------------------
103   !!   'key_dimgout'                                      DIMG output file
104   !!----------------------------------------------------------------------
105#   include "diawri_dimg.h90"
106
107#else
108   !!----------------------------------------------------------------------
109   !!   Default option                                   NetCDF output file
110   !!----------------------------------------------------------------------
111# if defined key_iomput
112   !!----------------------------------------------------------------------
113   !!   'key_iomput'                                        use IOM library
114   !!----------------------------------------------------------------------
115
116   SUBROUTINE dia_wri( kt )
117      !!---------------------------------------------------------------------
118      !!                  ***  ROUTINE dia_wri  ***
119      !!                   
120      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
121      !!      NETCDF format is used by default
122      !!
123      !! ** Method  :  use iom_put
124      !!----------------------------------------------------------------------
125      !!
126      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
127      !!
128      INTEGER                      ::   ji, jj, jk              ! dummy loop indices
129      INTEGER                      ::   jkbot                   !
130      REAL(wp)                     ::   zztmp, zztmpx, zztmpy, zmdi   !
131      !!
132      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
133      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
134      !!----------------------------------------------------------------------
135      !
136      zmdi=1.e+20 !missing data indicator for masking
137
138      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
139      !
140      CALL wrk_alloc( jpi , jpj      , z2d )
141      CALL wrk_alloc( jpi , jpj, jpk , z3d )
142      !
143      ! Output the initial state and forcings
144      IF( ninist == 1 ) THEN                       
145         CALL dia_wri_state( 'output.init', kt )
146         ninist = 0
147      ENDIF
148
149      IF( .NOT.lk_vvl ) THEN
150         CALL iom_put( "e3t" , fse3t_n(:,:,:) )
151         CALL iom_put( "e3u" , fse3u_n(:,:,:) )
152         CALL iom_put( "e3v" , fse3v_n(:,:,:) )
153         CALL iom_put( "e3w" , fse3w_n(:,:,:) )
154      ENDIF
155
156      z2d = sshn(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
157      CALL iom_put( "ssh" , z2d )                 ! sea surface height
158      if( iom_use('ssh2') )   CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
159     
160      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
161      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
162      IF ( iom_use("sbt") ) THEN
163         DO jj = 1, jpj
164            DO ji = 1, jpi
165               jkbot = mbkt(ji,jj)
166               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
167            END DO
168         END DO
169         CALL iom_put( "sbt", z2d )                ! bottom temperature
170      ENDIF
171     
172      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
173      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
174      IF ( iom_use("sbs") ) THEN
175         DO jj = 1, jpj
176            DO ji = 1, jpi
177               jkbot = mbkt(ji,jj)
178               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal)
179            END DO
180         END DO
181         CALL iom_put( "sbs", z2d )                ! bottom salinity
182      ENDIF
183
184      IF ( iom_use("taubot") ) THEN                ! bottom stress
185         z2d(:,:) = 0._wp
186         DO jj = 2, jpjm1
187            DO ji = fs_2, fs_jpim1   ! vector opt.
188               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  &
189                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )     
190               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
191                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
192               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
193               !
194            ENDDO
195         ENDDO
196         CALL lbc_lnk( z2d, 'T', 1. )
197         CALL iom_put( "taubot", z2d )           
198      ENDIF
199         
200      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current
201      z2d = un(:,:,1)*umask(:,:,1) + zmdi*(1.0-umask(:,:,1))
202      CALL iom_put(  "ssu", z2d         )    ! surface i-current
203      IF ( iom_use("sbu") ) THEN
204         DO jj = 1, jpj
205            DO ji = 1, jpi
206               jkbot = mbku(ji,jj)
207               z2d(ji,jj) = un(ji,jj,jkbot)
208            END DO
209         END DO
210         CALL iom_put( "sbu", z2d )                ! bottom i-current
211      ENDIF
212#if defined key_dynspg_ts
213      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current
214#else
215      CALL iom_put(  "ubar", un_b(:,:)        )    ! barotropic i-current
216#endif
217     
218      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
219      z2d = vn(:,:,1)*vmask(:,:,1) + zmdi*(1.0-vmask(:,:,1))     
220      CALL iom_put(  "ssv", z2d         )    ! surface j-current
221      IF ( iom_use("sbv") ) THEN
222         DO jj = 1, jpj
223            DO ji = 1, jpi
224               jkbot = mbkv(ji,jj)
225               z2d(ji,jj) = vn(ji,jj,jkbot)
226            END DO
227         END DO
228         CALL iom_put( "sbv", z2d )                ! bottom j-current
229      ENDIF
230#if defined key_dynspg_ts
231      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic j-current
232#else
233      CALL iom_put(  "vbar", vn_b(:,:)        )    ! barotropic j-current
234#endif
235
236      CALL iom_put( "woce", wn )                   ! vertical velocity
237      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
238         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
239         z2d(:,:) = rau0 * e12t(:,:)
240         DO jk = 1, jpk
241            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
242         END DO
243         CALL iom_put( "w_masstr" , z3d ) 
244         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
245      ENDIF
246
247      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef.
248      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef.
249      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm)
250
251      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
252         DO jj = 2, jpjm1                                    ! sst gradient
253            DO ji = fs_2, fs_jpim1   ! vector opt.
254               zztmp      = tsn(ji,jj,1,jp_tem)
255               zztmpx     = ( tsn(ji+1,jj  ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) / e1u(ji-1,jj  )
256               zztmpy     = ( tsn(ji  ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) / e2v(ji  ,jj-1)
257               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
258                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
259            END DO
260         END DO
261         CALL lbc_lnk( z2d, 'T', 1. )
262         CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient
263         z2d(:,:) = SQRT( z2d(:,:) )
264         CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient
265      ENDIF
266         
267      ! clem: heat and salt content
268      IF( iom_use("heatc") ) THEN
269         z2d(:,:)  = 0._wp 
270         DO jk = 1, jpkm1
271            DO jj = 1, jpj
272               DO ji = 1, jpi
273                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
274               END DO
275            END DO
276         END DO
277         CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2)
278      ENDIF
279
280      IF( iom_use("saltc") ) THEN
281         z2d(:,:)  = 0._wp 
282         DO jk = 1, jpkm1
283            DO jj = 1, jpj
284               DO ji = 1, jpi
285                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
286               END DO
287            END DO
288         END DO
289         CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2)
290      ENDIF
291      !
292      IF ( iom_use("eken") ) THEN
293         rke(:,:,jk) = 0._wp                               !      kinetic energy
294         DO jk = 1, jpkm1
295            DO jj = 2, jpjm1
296               DO ji = fs_2, fs_jpim1   ! vector opt.
297                  zztmp   = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
298                  zztmpx  = 0.5 * (  un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk)    &
299                     &             + un(ji  ,jj,jk) * un(ji  ,jj,jk) * e2u(ji  ,jj) * fse3u(ji  ,jj,jk) )  &
300                     &          *  zztmp 
301                  !
302                  zztmpy  = 0.5 * (  vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk)    &
303                     &             + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1v(ji,jj  ) * fse3v(ji,jj  ,jk) )  &
304                     &          *  zztmp 
305                  !
306                  rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy )
307                  !
308               ENDDO
309            ENDDO
310         ENDDO
311         CALL lbc_lnk( rke, 'T', 1. )
312         CALL iom_put( "eken", rke )           
313      ENDIF
314         
315      IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
316         z3d(:,:,jpk) = 0.e0
317         DO jk = 1, jpkm1
318            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
319         END DO
320         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
321      ENDIF
322     
323      IF( iom_use("u_heattr") ) THEN
324         z2d(:,:) = 0.e0 
325         DO jk = 1, jpkm1
326            DO jj = 2, jpjm1
327               DO ji = fs_2, fs_jpim1   ! vector opt.
328                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
329               END DO
330            END DO
331         END DO
332         CALL lbc_lnk( z2d, 'U', -1. )
333         CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction
334      ENDIF
335
336      IF( iom_use("u_salttr") ) THEN
337         z2d(:,:) = 0.e0 
338         DO jk = 1, jpkm1
339            DO jj = 2, jpjm1
340               DO ji = fs_2, fs_jpim1   ! vector opt.
341                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
342               END DO
343            END DO
344         END DO
345         CALL lbc_lnk( z2d, 'U', -1. )
346         CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction
347      ENDIF
348
349     
350      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
351         z3d(:,:,jpk) = 0.e0
352         DO jk = 1, jpkm1
353            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
354         END DO
355         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
356      ENDIF
357     
358      IF( iom_use("v_heattr") ) THEN
359         z2d(:,:) = 0.e0 
360         DO jk = 1, jpkm1
361            DO jj = 2, jpjm1
362               DO ji = fs_2, fs_jpim1   ! vector opt.
363                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
364               END DO
365            END DO
366         END DO
367         CALL lbc_lnk( z2d, 'V', -1. )
368         CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction
369      ENDIF
370
371      IF( iom_use("v_salttr") ) THEN
372         z2d(:,:) = 0.e0 
373         DO jk = 1, jpkm1
374            DO jj = 2, jpjm1
375               DO ji = fs_2, fs_jpim1   ! vector opt.
376                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
377               END DO
378            END DO
379         END DO
380         CALL lbc_lnk( z2d, 'V', -1. )
381         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction
382      ENDIF
383      !
384      CALL wrk_dealloc( jpi , jpj      , z2d )
385      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
386      !
387      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
388      !
389   END SUBROUTINE dia_wri
390
391#else
392   !!----------------------------------------------------------------------
393   !!   Default option                                  use IOIPSL  library
394   !!----------------------------------------------------------------------
395
396   SUBROUTINE dia_wri( kt )
397      !!---------------------------------------------------------------------
398      !!                  ***  ROUTINE dia_wri  ***
399      !!                   
400      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
401      !!      NETCDF format is used by default
402      !!
403      !! ** Method  :   At the beginning of the first time step (nit000),
404      !!      define all the NETCDF files and fields
405      !!      At each time step call histdef to compute the mean if ncessary
406      !!      Each nwrite time step, output the instantaneous or mean fields
407      !!----------------------------------------------------------------------
408      !!
409      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
410      !!
411      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
412      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
413      INTEGER  ::   inum = 11                                ! temporary logical unit
414      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
415      INTEGER  ::   ierr                                     ! error code return from allocation
416      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
417      INTEGER  ::   jn, ierror                               ! local integers
418      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt           ! local scalars
419      !!
420      REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
421      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
422      !!----------------------------------------------------------------------
423      !
424      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
425      !
426      CALL wrk_alloc( jpi , jpj      , zw2d )
427      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_alloc( jpi , jpj , jpk  , zw3d )
428      !
429      ! Output the initial state and forcings
430      IF( ninist == 1 ) THEN                       
431         CALL dia_wri_state( 'output.init', kt )
432         ninist = 0
433      ENDIF
434      !
435      ! 0. Initialisation
436      ! -----------------
437
438      ! local variable for debugging
439      ll_print = .FALSE.
440      ll_print = ll_print .AND. lwp
441
442      ! Define frequency of output and means
443      zdt = rdt
444      IF( nacc == 1 ) zdt = rdtmin
445      clop = "x"         ! no use of the mask value (require less cpu time, and otherwise the model crashes)
446#if defined key_diainstant
447      zsto = nwrite * zdt
448      clop = "inst("//TRIM(clop)//")"
449#else
450      zsto=zdt
451      clop = "ave("//TRIM(clop)//")"
452#endif
453      zout = nwrite * zdt
454      zmax = ( nitend - nit000 + 1 ) * zdt
455
456      ! Define indices of the horizontal output zoom and vertical limit storage
457      iimi = 1      ;      iima = jpi
458      ijmi = 1      ;      ijma = jpj
459      ipk = jpk
460
461      ! define time axis
462      it = kt
463      itmod = kt - nit000 + 1
464
465
466      ! 1. Define NETCDF files and fields at beginning of first time step
467      ! -----------------------------------------------------------------
468
469      IF( kt == nit000 ) THEN
470
471         ! Define the NETCDF files (one per grid)
472
473         ! Compute julian date from starting date of the run
474         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
475         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
476         IF(lwp)WRITE(numout,*)
477         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
478            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
479         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
480                                 ' limit storage in depth = ', ipk
481
482         ! WRITE root name in date.file for use by postpro
483         IF(lwp) THEN
484            CALL dia_nam( clhstnam, nwrite,' ' )
485            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
486            WRITE(inum,*) clhstnam
487            CLOSE(inum)
488         ENDIF
489
490         ! Define the T grid FILE ( nid_T )
491
492         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
493         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
494         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
495            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
496            &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
497         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
498            &           "m", ipk, gdept_1d, nz_T, "down" )
499         !                                                            ! Index of ocean points
500         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
501         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
502         !
503         IF( ln_icebergs ) THEN
504            !
505            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
506            !! that routine is called from nemogcm, so do it here immediately before its needed
507            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
508            IF( lk_mpp )   CALL mpp_sum( ierror )
509            IF( ierror /= 0 ) THEN
510               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
511               RETURN
512            ENDIF
513            !
514            !! iceberg vertical coordinate is class number
515            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
516               &           "number", nclasses, class_num, nb_T )
517            !
518            !! each class just needs the surface index pattern
519            ndim_bT = 3
520            DO jn = 1,nclasses
521               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
522            ENDDO
523            !
524         ENDIF
525
526         ! Define the U grid FILE ( nid_U )
527
528         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
529         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
530         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
531            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
532            &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
533         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
534            &           "m", ipk, gdept_1d, nz_U, "down" )
535         !                                                            ! Index of ocean points
536         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
537         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
538
539         ! Define the V grid FILE ( nid_V )
540
541         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
542         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
543         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
544            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
545            &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
546         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
547            &          "m", ipk, gdept_1d, nz_V, "down" )
548         !                                                            ! Index of ocean points
549         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
550         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
551
552         ! Define the W grid FILE ( nid_W )
553
554         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
555         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
556         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
557            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
558            &          nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
559         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
560            &          "m", ipk, gdepw_1d, nz_W, "down" )
561
562
563         ! Declare all the output fields as NETCDF variables
564
565         !                                                                                      !!! nid_T : 3D
566         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
567            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
568         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
569            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
570         IF(  lk_vvl  ) THEN
571            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
572            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
573            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
574            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
575            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
576            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
577         ENDIF
578         !                                                                                      !!! nid_T : 2D
579         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
580            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
581         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
582            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
583         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
584            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
585         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
586            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
587         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
588            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
589         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
590            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
591         IF(  .NOT. lk_vvl  ) THEN
592            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
593            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
594            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
595            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
596            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
597            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
598         ENDIF
599         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
600            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
601         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
602            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
603         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
604            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
605         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
606            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
607         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
608            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
609         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
610            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
611!
612         IF( ln_icebergs ) THEN
613            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
614               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
615            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
616               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
617            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
618               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
619            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
620               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
621            IF( ln_bergdia ) THEN
622               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
623                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
624               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
625                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
626               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
627                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
628               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
629                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
630               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
631                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
632               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
633                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
634               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
635                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
636               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
637                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
638               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
639                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
640               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
641                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
642            ENDIF
643         ENDIF
644
645         IF( .NOT. ln_cpl ) THEN
646            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
647               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
648            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
649               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
650            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
651               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
652         ENDIF
653
654         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
655            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
656               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
657            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
658               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
659            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
660               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
661         ENDIF
662         
663         clmx ="l_max(only(x))"    ! max index on a period
664         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
665            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
666#if defined key_diahth
667         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
668            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
669         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
670            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
671         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
672            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
673         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "W"      ,   & ! htc3
674            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
675#endif
676
677         IF( ln_cpl .AND. nn_ice == 2 ) THEN
678            CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
679               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
680            CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
681               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
682         ENDIF
683
684         CALL histend( nid_T, snc4chunks=snc4set )
685
686         !                                                                                      !!! nid_U : 3D
687         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
688            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
689         IF( ln_traldf_gdia ) THEN
690            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
691                 &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
692         ELSE
693#if defined key_diaeiv
694            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
695            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
696#endif
697         END IF
698         !                                                                                      !!! nid_U : 2D
699         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
700            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
701
702         CALL histend( nid_U, snc4chunks=snc4set )
703
704         !                                                                                      !!! nid_V : 3D
705         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
706            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
707         IF( ln_traldf_gdia ) THEN
708            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
709                 &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
710         ELSE 
711#if defined key_diaeiv
712            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
713            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
714#endif
715         END IF
716         !                                                                                      !!! nid_V : 2D
717         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
718            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
719
720         CALL histend( nid_V, snc4chunks=snc4set )
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         IF( ln_traldf_gdia ) THEN
726            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
727                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
728         ELSE
729#if defined key_diaeiv
730            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
731                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
732#endif
733         END IF
734         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
735            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
736         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
737            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
738
739         IF( lk_zdfddm ) THEN
740            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
741               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
742         ENDIF
743         !                                                                                      !!! nid_W : 2D
744#if defined key_traldf_c2d
745         CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity"           , "m2/s"   ,   &  ! ahtw
746            &          jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
747# if defined key_traldf_eiv 
748            CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s",   &  ! aeiw
749               &       jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
750# endif
751#endif
752
753         CALL histend( nid_W, snc4chunks=snc4set )
754
755         IF(lwp) WRITE(numout,*)
756         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
757         IF(ll_print) CALL FLUSH(numout )
758
759      ENDIF
760
761      ! 2. Start writing data
762      ! ---------------------
763
764      ! ndex(1) est utilise ssi l'avant dernier argument est different de
765      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
766      ! donne le nombre d'elements, et ndex la liste des indices a sortir
767
768      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
769         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
770         WRITE(numout,*) '~~~~~~ '
771      ENDIF
772
773      IF( lk_vvl ) THEN
774         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
775         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
776         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
777         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
778      ELSE
779         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
780         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
781         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
782         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
783      ENDIF
784      IF( lk_vvl ) THEN
785         zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
786         CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
787         CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
788         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
789      ENDIF
790      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
791      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
792      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
793      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
794                                                                                  ! (includes virtual salt flux beneath ice
795                                                                                  ! in linear free surface case)
796      IF( .NOT. lk_vvl ) THEN
797         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
798         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
799         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
800         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
801      ENDIF
802      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
803      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
804      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
805      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
806      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
807      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
808!
809      IF( ln_icebergs ) THEN
810         !
811         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
812         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
813         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
814         !
815         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
816         !
817         IF( ln_bergdia ) THEN
818            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
819            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
820            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
821            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
822            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
823            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
824            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
825            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
826            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
827            !
828            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
829         ENDIF
830      ENDIF
831
832      IF( .NOT. ln_cpl ) THEN
833         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
834         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
835         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
836         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
837      ENDIF
838      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
839         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
840         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
841         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
842         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
843      ENDIF
844!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
845!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
846
847#if defined key_diahth
848      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
849      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
850      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
851      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
852#endif
853
854      IF( ln_cpl .AND. nn_ice == 2 ) THEN
855         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
856         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
857      ENDIF
858
859      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
860      IF( ln_traldf_gdia ) THEN
861         IF (.not. ALLOCATED(psix_eiv))THEN
862            ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
863            IF( lk_mpp   )   CALL mpp_sum ( ierr )
864            IF( ierr > 0 )   CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
865            psix_eiv(:,:,:) = 0.0_wp
866            psiy_eiv(:,:,:) = 0.0_wp
867         ENDIF
868         DO jk=1,jpkm1
869            zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz
870         END DO
871         zw3d(:,:,jpk) = 0._wp
872         CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U )           ! i-eiv current
873      ELSE
874#if defined key_diaeiv
875         CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U )          ! i-eiv current
876#endif
877      ENDIF
878      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
879
880      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
881      IF( ln_traldf_gdia ) THEN
882         DO jk=1,jpk-1
883            zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz
884         END DO
885         zw3d(:,:,jpk) = 0._wp
886         CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V )           ! j-eiv current
887      ELSE
888#if defined key_diaeiv
889         CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V )          ! j-eiv current
890#endif
891      ENDIF
892      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
893
894      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
895      IF( ln_traldf_gdia ) THEN
896         DO jk=1,jpk-1
897            DO jj = 2, jpjm1
898               DO ji = fs_2, fs_jpim1  ! vector opt.
899                  zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
900                       &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
901               END DO
902            END DO
903         END DO
904         zw3d(:,:,jpk) = 0._wp
905         CALL histwrite( nid_W, "voveeivw", it, zw3d          , ndim_T, ndex_T )    ! vert. eiv current
906      ELSE
907#   if defined key_diaeiv
908         CALL histwrite( nid_W, "voveeivw", it, w_eiv          , ndim_T, ndex_T )    ! vert. eiv current
909#   endif
910      ENDIF
911      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
912      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
913      IF( lk_zdfddm ) THEN
914         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
915      ENDIF
916#if defined key_traldf_c2d
917      CALL histwrite( nid_W, "soleahtw", it, ahtw          , ndim_hT, ndex_hT )   ! lateral eddy diff. coef.
918# if defined key_traldf_eiv
919         CALL histwrite( nid_W, "soleaeiw", it, aeiw       , ndim_hT, ndex_hT )   ! EIV coefficient at w-point
920# endif
921#endif
922
923      ! 3. Close all files
924      ! ---------------------------------------
925      IF( kt == nitend ) THEN
926         CALL histclo( nid_T )
927         CALL histclo( nid_U )
928         CALL histclo( nid_V )
929         CALL histclo( nid_W )
930      ENDIF
931      !
932      CALL wrk_dealloc( jpi , jpj      , zw2d )
933      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_dealloc( jpi , jpj , jpk  , zw3d )
934      !
935      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
936      !
937   END SUBROUTINE dia_wri
938# endif
939
940#endif
941
942   SUBROUTINE dia_wri_state( cdfile_name, kt )
943      !!---------------------------------------------------------------------
944      !!                 ***  ROUTINE dia_wri_state  ***
945      !!       
946      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
947      !!      the instantaneous ocean state and forcing fields.
948      !!        Used to find errors in the initial state or save the last
949      !!      ocean state in case of abnormal end of a simulation
950      !!
951      !! ** Method  :   NetCDF files using ioipsl
952      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
953      !!      File 'output.abort.nc' is created in case of abnormal job end
954      !!----------------------------------------------------------------------
955      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
956      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
957      !!
958      CHARACTER (len=32) :: clname
959      CHARACTER (len=40) :: clop
960      INTEGER  ::   id_i , nz_i, nh_i       
961      INTEGER, DIMENSION(1) ::   idex             ! local workspace
962      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt
963      !!----------------------------------------------------------------------
964      !
965!     IF( nn_timing == 1 )   CALL timing_start('dia_wri_state') ! not sure this works for routines not called in first timestep
966
967      ! 0. Initialisation
968      ! -----------------
969
970      ! Define name, frequency of output and means
971      clname = cdfile_name
972      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
973      zdt  = rdt
974      zsto = rdt
975      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
976      zout = rdt
977      zmax = ( nitend - nit000 + 1 ) * zdt
978
979      IF(lwp) WRITE(numout,*)
980      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
981      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
982      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
983
984
985      ! 1. Define NETCDF files and fields at beginning of first time step
986      ! -----------------------------------------------------------------
987
988      ! Compute julian date from starting date of the run
989      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
990      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
991      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
992          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
993      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
994          "m", jpk, gdept_1d, nz_i, "down")
995
996      ! Declare all the output fields as NetCDF variables
997
998      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
999         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1000      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
1001         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1002      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
1003         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
1004      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
1005         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1006      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
1007         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1008      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
1009         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1010      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
1011         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1012      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
1013         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1014      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
1015         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1016      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
1017         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1018      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
1019         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1020      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
1021         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1022      IF( lk_vvl ) THEN
1023         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      ,   &   ! t-point depth
1024            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1025      END IF
1026
1027#if defined key_lim2
1028      CALL lim_wri_state_2( kt, id_i, nh_i )
1029#elif defined key_lim3
1030      CALL lim_wri_state( kt, id_i, nh_i )
1031#else
1032      CALL histend( id_i, snc4chunks=snc4set )
1033#endif
1034
1035      ! 2. Start writing data
1036      ! ---------------------
1037      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
1038      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
1039      ! donne le nombre d'elements, et idex la liste des indices a sortir
1040      idex(1) = 1   ! init to avoid compil warning
1041
1042      ! Write all fields on T grid
1043      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
1044      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
1045      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
1046      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
1047      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
1048      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
1049      CALL histwrite( id_i, "sowaflup", kt, (emp-rnf )       , jpi*jpj    , idex )    ! freshwater budget
1050      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
1051      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
1052      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
1053      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
1054      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
1055
1056      ! 3. Close the file
1057      ! -----------------
1058      CALL histclo( id_i )
1059#if ! defined key_iomput && ! defined key_dimgout
1060      IF( ninist /= 1  ) THEN
1061         CALL histclo( nid_T )
1062         CALL histclo( nid_U )
1063         CALL histclo( nid_V )
1064         CALL histclo( nid_W )
1065      ENDIF
1066#endif
1067       
1068!     IF( nn_timing == 1 )   CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep
1069      !
1070
1071   END SUBROUTINE dia_wri_state
1072   !!======================================================================
1073END MODULE diawri
Note: See TracBrowser for help on using the repository browser.