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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4990

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

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

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