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

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6010

Last change on this file since 6010 was 6010, checked in by timgraham, 9 years ago

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