New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5239

Last change on this file since 5239 was 5239, checked in by davestorkey, 9 years ago

Commit changes to UKMO nn_etau revision branch (including clearing
SVN keywords).

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