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/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5676

Last change on this file since 5676 was 5676, checked in by jwhile, 9 years ago

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