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 on Ticket #465 – Attachment – NEMO

Ticket #465: diawri.F90

File diawri.F90, 35.0 KB (added by rachel.furner, 14 years ago)
Line 
1MODULE diawri
2   !!======================================================================
3   !!                     ***  MODULE  diawri  ***
4   !! Ocean diagnostics :  write ocean output files
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !! * Modules used
9   USE oce             ! ocean dynamics and tracers
10   USE dom_oce         ! ocean space and time domain
11   USE zdf_oce         ! ocean vertical physics
12   USE ldftra_oce      ! ocean active tracers: lateral physics
13   USE ldfdyn_oce      ! ocean dynamics: lateral physics
14   USE sol_oce         ! solver variables
15   USE sbc_oce         ! Surface boundary condition: ocean fields
16   USE sbc_ice         ! Surface boundary condition: ice fields
17   USE sbcssr          ! restoring term toward SST/SSS climatology
18   USE phycst          ! physical constants
19   USE zdfmxl          ! mixed layer
20   USE daymod          ! calendar
21   USE dianam          ! build name of file (routine)
22   USE zdfddm          ! vertical  physics: double diffusion
23   USE diahth          ! thermocline diagnostics
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE in_out_manager  ! I/O manager
26   USE diadimg         ! dimg direct access file format output
27   USE iom
28   USE ioipsl
29#if defined key_lim2
30   USE limwri_2 
31#endif
32   IMPLICIT NONE
33   PRIVATE
34
35   !! * Accessibility
36   PUBLIC dia_wri                 ! routines called by step.F90
37   PUBLIC dia_wri_state
38
39   !! * Module variables
40   INTEGER ::   &
41      nid_T, nz_T, nh_T, ndim_T, ndim_hT,      &   ! grid_T file
42      nid_U, nz_U, nh_U, ndim_U, ndim_hU,      &   ! grid_U file
43      nid_V, nz_V, nh_V, ndim_V, ndim_hV,      &   ! grid_V file
44      nid_W, nz_W, nh_W,                       &   ! grid_W file
45      ndex(1)                                      ! ???
46   INTEGER, DIMENSION(jpi*jpj) ::   &
47      ndex_hT, ndex_hU, ndex_hV
48   INTEGER, DIMENSION(jpi*jpj*jpk) ::   &
49      ndex_T, ndex_U, ndex_V
50
51   !! * Substitutions
52#  include "zdfddm_substitute.h90"
53   !!----------------------------------------------------------------------
54   !!   OPA 9.0 , LOCEAN-IPSL (2005)
55   !! $Id: diawri.F90 1585 2009-08-06 07:29:43Z smasson $
56   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
57   !!----------------------------------------------------------------------
58
59CONTAINS
60
61#if defined key_dimgout
62   !!----------------------------------------------------------------------
63   !!   dia_wri       : create the dimg direct access output file (mpp)
64   !!----------------------------------------------------------------------
65#   include "diawri_dimg.h90"
66
67#else
68   !!----------------------------------------------------------------------
69   !!   Default option                                   NetCDF output file
70   !!----------------------------------------------------------------------
71   !!   dia_wri       : create the standart NetCDF output files
72   !!   dia_wri_state : create an output NetCDF file for a single
73   !!                   instantaeous ocean state and forcing fields
74   !!----------------------------------------------------------------------
75# if defined key_iomput
76   SUBROUTINE dia_wri( kt )
77      !!---------------------------------------------------------------------
78      !!                  ***  ROUTINE dia_wri  ***
79      !!                   
80      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
81      !!      NETCDF format is used by default
82      !!
83      !! ** Method  :  use iom_put
84      !!
85      !! History :
86      !!   3.2  !  05-11  (B. Lemaire) creation from old diawri
87      !!----------------------------------------------------------------------
88      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
89      !!----------------------------------------------------------------------
90      !
91      ! Output the initial state and forcings
92      IF( ninist == 1 ) THEN                       
93         CALL dia_wri_state( 'output.init', kt )
94         ninist = 0
95      ENDIF
96
97      CALL iom_put( "toce"   , tn        )    ! temperature
98      CALL iom_put( "soce"   , sn        )    ! salinity
99      CALL iom_put( "sst"    , tn(:,:,1) )    ! sea surface temperature
100      CALL iom_put( "sss"    , sn(:,:,1) )    ! sea surface salinity
101      CALL iom_put( "uoce"   , un        )    ! i-current     
102      CALL iom_put( "voce"   , vn        )    ! j-current
103     
104      CALL iom_put( "avt"    , avt       )    ! T vert. eddy diff. coef.
105      CALL iom_put( "avm"    , avmu      )    ! T vert. eddy visc. coef.
106      IF( lk_zdfddm ) THEN
107         CALL iom_put( "avs", fsavs(:,:,:) )    ! S vert. eddy diff. coef.
108      ENDIF
109
110   END SUBROUTINE dia_wri
111
112#else
113   SUBROUTINE dia_wri( kt )
114      !!---------------------------------------------------------------------
115      !!                  ***  ROUTINE dia_wri  ***
116      !!                   
117      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
118      !!      NETCDF format is used by default
119      !!
120      !! ** Method  :   At the beginning of the first time step (nit000),
121      !!      define all the NETCDF files and fields
122      !!      At each time step call histdef to compute the mean if ncessary
123      !!      Each nwrite time step, output the instantaneous or mean fields
124      !!
125      !! History :
126      !!        !  91-03  (M.-A. Foujols)  Original code
127      !!        !  91-11  (G. Madec)
128      !!        !  92-06  (M. Imbard)  correction restart file
129      !!        !  92-07  (M. Imbard)  split into diawri and rstwri
130      !!        !  93-03  (M. Imbard)  suppress writibm
131      !!        !  98-01  (C. Levy)  NETCDF format using ioipsl INTERFACE
132      !!        !  99-02  (E. Guilyardi)  name of netCDF files + variables
133      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
134      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
135      !!----------------------------------------------------------------------
136      !! * Arguments
137      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
138
139      !! * Local declarations
140      LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout
141      CHARACTER (len=40) ::           &
142         clhstnam, clop, clmx            ! temporary names
143      INTEGER ::   inum = 11             ! temporary logical unit
144      INTEGER ::   &
145         iimi, iima, ipk, it, itmod,  &  ! temporary integers
146         ijmi, ijma                      !    "          "
147      REAL(wp) ::   &
148         zsto, zout, zmax,            &  ! temporary scalars
149         zjulian, zdt                    !    "         "
150      REAL(wp), DIMENSION(jpi,jpj) :: &
151         zw2d                            ! temporary workspace
152      !!----------------------------------------------------------------------
153      !
154      ! Output the initial state and forcings
155      IF( ninist == 1 ) THEN                       
156         CALL dia_wri_state( 'output.init', kt )
157         ninist = 0
158      ENDIF
159      !
160      ! 0. Initialisation
161      ! -----------------
162
163      ! local variable for debugging
164      ll_print = .FALSE.
165      ll_print = ll_print .AND. lwp
166
167      ! Define frequency of output and means
168      zdt = rdt
169      IF( nacc == 1 ) zdt = rdtmin
170      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
171      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
172      ENDIF
173#if defined key_diainstant
174      zsto = nwrite * zdt
175      clop = "inst("//TRIM(clop)//")"
176#else
177      zsto=zdt
178      clop = "ave("//TRIM(clop)//")"
179#endif
180      zout = nwrite * zdt
181      zmax = ( nitend - nit000 + 1 ) * zdt
182
183      ! Define indices of the horizontal output zoom and vertical limit storage
184      iimi = 1      ;      iima = jpi
185      ijmi = 1      ;      ijma = jpj
186      ipk = jpk
187
188      ! define time axis
189      it = kt
190      itmod = kt - nit000 + 1
191
192
193      ! 1. Define NETCDF files and fields at beginning of first time step
194      ! -----------------------------------------------------------------
195
196      IF( kt == nit000 ) THEN
197
198         ! Define the NETCDF files (one per grid)
199
200         ! Compute julian date from starting date of the run
201         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
202         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
203         IF(lwp)WRITE(numout,*)
204         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
205            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
206         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
207                                 ' limit storage in depth = ', ipk
208
209         ! WRITE root name in date.file for use by postpro
210         IF(lwp) THEN
211            CALL dia_nam( clhstnam, nwrite,' ' )
212            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
213            WRITE(inum,*) clhstnam
214            CLOSE(inum)
215         ENDIF
216
217         ! Define the T grid FILE ( nid_T )
218
219         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
220         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
221         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
222            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
223            &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom )
224         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
225            &           "m", ipk, gdept_0, nz_T, "down" )
226         !                                                            ! Index of ocean points
227         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
228         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
229
230         ! Define the U grid FILE ( nid_U )
231
232         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
233         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
234         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
235            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
236            &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom )
237         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
238            &           "m", ipk, gdept_0, nz_U, "down" )
239         !                                                            ! Index of ocean points
240         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
241         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
242
243         ! Define the V grid FILE ( nid_V )
244
245         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
246         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
247         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
248            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
249            &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom )
250         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
251            &          "m", ipk, gdept_0, nz_V, "down" )
252         !                                                            ! Index of ocean points
253         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
254         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
255
256         ! Define the W grid FILE ( nid_W )
257
258         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
259         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
260         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
261            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
262            &          nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom )
263         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
264            &          "m", ipk, gdepw_0, nz_W, "down" )
265
266
267         ! Declare all the output fields as NETCDF variables
268
269         !                                                                                      !!! nid_T : 3D
270         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
271            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
272         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
273            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
274         !                                                                                      !!! nid_T : 2D
275         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
276            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
277         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
278            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
279         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
280            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
281!!$#if defined key_lim3 || defined key_lim2
282!!$         ! sowaflup = sowaflep + sorunoff + sowafldp + a term associated to
283!!$         !    internal damping to Levitus that can be diagnosed from others
284!!$         ! sowaflcd = sowaflep + sorunoff + sowafldp + iowaflup
285!!$         CALL histdef( nid_T, "iowaflup", "Ice=>ocean net freshwater"          , "kg/m2/s",   &  ! fsalt
286!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
287!!$         CALL histdef( nid_T, "sowaflep", "atmos=>ocean net freshwater"        , "kg/m2/s",   &  ! fmass
288!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
289!!$#endif
290         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
291            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
292!!$         CALL histdef( nid_T, "sorunoff", "Runoffs"                            , "Kg/m2/s",   &  ! runoffs
293!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
294         CALL histdef( nid_T, "sowaflcd", "concentration/dilution water flux"  , "kg/m2/s",   &  ! (emps-rnf)
295            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
296         CALL histdef( nid_T, "sosalflx", "Surface Salt Flux"                  , "Kg/m2/s",   &  ! (emps-rnf) * sn
297            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
298         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
299            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
300         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
301            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
302         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
303            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
304         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
305            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
306         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
307            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
308#if ! defined key_coupled
309         CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
310            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
311         CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
312            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
313         CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
314            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
315#endif
316
317
318
319#if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 )
320         CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
321            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
322         CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
323            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
324         CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
325            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
326#endif
327#if defined key_diaspr
328         CALL histdef( nid_T, "sosurfps", "Surface Pressure"                   , "cm"     ,   &  ! sp
329            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
330#endif
331         clmx ="l_max(only(x))"    ! max index on a period
332         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
333            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
334#if defined key_diahth
335         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
336            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
337         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
338            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
339         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
340            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
341         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "W"      ,   & ! htc3
342            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
343#endif
344
345#if defined key_coupled 
346# if defined key_lim3
347         Must be adapted to LIM3
348# else
349         CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
350            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
351         CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
352            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
353# endif 
354#endif
355
356         CALL histend( nid_T )
357
358         !                                                                                      !!! nid_U : 3D
359         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
360            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
361#if defined key_diaeiv
362         CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
363            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
364#endif
365         !                                                                                      !!! nid_U : 2D
366         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
367            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
368
369         CALL histend( nid_U )
370
371         !                                                                                      !!! nid_V : 3D
372         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
373            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
374#if defined key_diaeiv
375         CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
376            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
377#endif
378         !                                                                                      !!! nid_V : 2D
379         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
380            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
381
382         CALL histend( nid_V )
383
384         !                                                                                      !!! nid_W : 3D
385         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
386            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
387#if defined key_diaeiv
388         CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
389            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
390#endif
391         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
392            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
393         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
394            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
395
396         IF( lk_zdfddm ) THEN
397            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
398               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
399         ENDIF
400         !                                                                                      !!! nid_W : 2D
401#if defined key_traldf_c2d
402         CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity"           , "m2/s"   ,   &  ! ahtw
403            &          jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
404# if defined key_traldf_eiv 
405            CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s",   &  ! aeiw
406               &       jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
407# endif
408#endif
409
410         CALL histend( nid_W )
411
412         IF(lwp) WRITE(numout,*)
413         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
414         IF(ll_print) CALL FLUSH(numout )
415
416      ENDIF
417
418      ! 2. Start writing data
419      ! ---------------------
420
421      ! ndex(1) est utilise ssi l'avant dernier argument est diffferent de
422      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
423      ! donne le nombre d'elements, et ndex la liste des indices a sortir
424
425      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
426         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
427         WRITE(numout,*) '~~~~~~ '
428      ENDIF
429
430      ! Write fields on T grid
431      CALL histwrite( nid_T, "votemper", it, tn            , ndim_T , ndex_T  )   ! temperature
432      CALL histwrite( nid_T, "vosaline", it, sn            , ndim_T , ndex_T  )   ! salinity
433      CALL histwrite( nid_T, "sosstsst", it, tn(:,:,1)     , ndim_hT, ndex_hT )   ! sea surface temperature
434      CALL histwrite( nid_T, "sosaline", it, sn(:,:,1)     , ndim_hT, ndex_hT )   ! sea surface salinity
435      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
436!!$#if  defined key_lim3 || defined key_lim2
437!!$      CALL histwrite( nid_T, "iowaflup", it, fsalt(:,:)    , ndim_hT, ndex_hT )   ! ice=>ocean water flux
438!!$      CALL histwrite( nid_T, "sowaflep", it, fmass(:,:)    , ndim_hT, ndex_hT )   ! atmos=>ocean water flux
439!!$#endif
440      CALL histwrite( nid_T, "sowaflup", it, ( emp - rnf )    , ndim_hT, ndex_hT )   ! upward water flux
441!!$      CALL histwrite( nid_T, "sorunoff", it, runoff        , ndim_hT, ndex_hT )   ! runoff
442      CALL histwrite( nid_T, "sowaflcd", it, ( emps - rnf )   , ndim_hT, ndex_hT )   ! c/d water flux
443      zw2d(:,:) = ( emps(:,:) - rnf(:,:) ) * sn(:,:,1) * tmask(:,:,1)
444      CALL histwrite( nid_T, "sosalflx", it, zw2d          , ndim_hT, ndex_hT )   ! c/d salt flux
445      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
446      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
447      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
448      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
449      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
450#if ! defined key_coupled
451      CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
452      CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
453      zw2d(:,:) = erp(:,:) * sn(:,:,1) * tmask(:,:,1)
454      CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
455#endif
456#if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 )
457      CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
458      CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
459         zw2d(:,:) = erp(:,:) * sn(:,:,1) * tmask(:,:,1)
460      CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
461#endif
462#if defined key_diaspr
463      CALL histwrite( nid_T, "sosurfps", it, gps           , ndim_hT, ndex_hT )   ! surface pressure
464#endif
465         zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
466      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
467
468#if defined key_diahth
469      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
470      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
471      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
472      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
473#endif
474
475#if defined key_coupled 
476# if defined key_lim3
477      Must be adapted for LIM3
478      CALL histwrite( nid_T, "soicetem", it, tn_ice        , ndim_hT, ndex_hT )   ! surf. ice temperature
479      CALL histwrite( nid_T, "soicealb", it, alb_ice       , ndim_hT, ndex_hT )   ! ice albedo
480# else
481      CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
482      CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
483# endif
484#endif
485         ! Write fields on U grid
486      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
487#if defined key_diaeiv
488      CALL histwrite( nid_U, "vozoeivu", it, u_eiv         , ndim_U , ndex_U )    ! i-eiv current
489#endif
490      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
491
492         ! Write fields on V grid
493      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
494#if defined key_diaeiv
495      CALL histwrite( nid_V, "vomeeivv", it, v_eiv         , ndim_V , ndex_V  )   ! j-eiv current
496#endif
497      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
498
499         ! Write fields on W grid
500      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
501#   if defined key_diaeiv
502      CALL histwrite( nid_W, "voveeivw", it, w_eiv          , ndim_T, ndex_T )    ! vert. eiv current
503#   endif
504      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
505      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
506      IF( lk_zdfddm ) THEN
507         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
508      ENDIF
509#if defined key_traldf_c2d
510      CALL histwrite( nid_W, "soleahtw", it, ahtw          , ndim_hT, ndex_hT )   ! lateral eddy diff. coef.
511# if defined key_traldf_eiv
512         CALL histwrite( nid_W, "soleaeiw", it, aeiw       , ndim_hT, ndex_hT )   ! EIV coefficient at w-point
513# endif
514#endif
515
516      ! 3. Close all files
517      ! ---------------------------------------
518      IF( kt == nitend ) THEN
519         CALL histclo( nid_T )
520         CALL histclo( nid_U )
521         CALL histclo( nid_V )
522         CALL histclo( nid_W )
523      ENDIF
524
525   END SUBROUTINE dia_wri
526# endif
527
528#endif
529
530   SUBROUTINE dia_wri_state( cdfile_name, kt )
531      !!---------------------------------------------------------------------
532      !!                 ***  ROUTINE dia_wri_state  ***
533      !!       
534      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
535      !!      the instantaneous ocean state and forcing fields.
536      !!        Used to find errors in the initial state or save the last
537      !!      ocean state in case of abnormal end of a simulation
538      !!
539      !! ** Method  :   NetCDF files using ioipsl
540      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
541      !!      File 'output.abort.nc' is created in case of abnormal job end
542      !!
543      !! History :
544      !!   8.2  !  00-06  (M. Imbard)  Original code (diabort.F)
545      !!   8.5  !  02-06  (A.Bozec, E. Durand)  Original code (diainit.F)
546      !!   9.0  !  02-12  (G. Madec)  merge of diabort and diainit, F90
547      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
548      !!----------------------------------------------------------------------
549      !! * Arguments
550      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
551      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
552
553      !! * Local declarations
554      CHARACTER (len=32) :: clname
555      CHARACTER (len=40) :: clop
556      INTEGER  ::   &
557         id_i , nz_i, nh_i       
558      INTEGER, DIMENSION(1) ::   &
559         idex             ! temprary workspace
560      REAL(wp) ::   &
561         zsto, zout, zmax,   &
562         zjulian, zdt
563      !!----------------------------------------------------------------------
564
565      ! 0. Initialisation
566      ! -----------------
567
568      ! Define name, frequency of output and means
569      clname = cdfile_name
570#if defined key_agrif
571      if ( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
572#endif
573      zdt  = rdt
574      zsto = rdt
575      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
576      zout = rdt
577      zmax = ( nitend - nit000 + 1 ) * zdt
578
579      IF(lwp) WRITE(numout,*)
580      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
581      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
582      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
583
584
585      ! 1. Define NETCDF files and fields at beginning of first time step
586      ! -----------------------------------------------------------------
587
588      ! Compute julian date from starting date of the run
589      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
590      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
591      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
592          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom )          ! Horizontal grid : glamt and gphit
593      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
594          "m", jpk, gdept_0, nz_i, "down")
595
596      ! Declare all the output fields as NetCDF variables
597
598      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
599         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
600      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
601         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
602      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
603         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
604      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
605         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
606      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
607         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
608      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
609         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
610      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
611         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
612      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
613         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
614      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
615         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
616      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
617         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
618      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
619         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
620      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
621         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
622
623#if defined key_lim2
624      CALL lim_wri_state_2( kt, id_i, nh_i )
625#else
626      CALL histend( id_i )
627#endif
628
629      ! 2. Start writing data
630      ! ---------------------
631      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
632      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
633      ! donne le nombre d'elements, et idex la liste des indices a sortir
634      idex(1) = 1   ! init to avoid compil warning
635
636      ! Write all fields on T grid
637      CALL histwrite( id_i, "votemper", kt, tn      , jpi*jpj*jpk, idex )    ! now temperature
638      CALL histwrite( id_i, "vosaline", kt, sn      , jpi*jpj*jpk, idex )    ! now salinity
639      CALL histwrite( id_i, "sossheig", kt, sshn     , jpi*jpj    , idex )    ! sea surface height
640      CALL histwrite( id_i, "vozocrtx", kt, un       , jpi*jpj*jpk, idex )    ! now i-velocity
641      CALL histwrite( id_i, "vomecrty", kt, vn       , jpi*jpj*jpk, idex )    ! now j-velocity
642      CALL histwrite( id_i, "vovecrtz", kt, wn       , jpi*jpj*jpk, idex )    ! now k-velocity
643      CALL histwrite( id_i, "sowaflup", kt, (emp-rnf), jpi*jpj    , idex )    ! freshwater budget
644      CALL histwrite( id_i, "sohefldo", kt, qsr + qns, jpi*jpj    , idex )    ! total heat flux
645      CALL histwrite( id_i, "soshfldo", kt, qsr      , jpi*jpj    , idex )    ! solar heat flux
646      CALL histwrite( id_i, "soicecov", kt, fr_i     , jpi*jpj    , idex )    ! ice fraction
647      CALL histwrite( id_i, "sozotaux", kt, utau     , jpi*jpj    , idex )    ! i-wind stress
648      CALL histwrite( id_i, "sometauy", kt, vtau     , jpi*jpj    , idex )    ! j-wind stress
649
650      ! 3. Close the file
651      ! -----------------
652      CALL histclo( id_i )
653#if ! defined key_iomput && ! defined key_dimgout
654      IF( ninist /= 1  ) THEN
655         CALL histclo( nid_T )
656         CALL histclo( nid_U )
657         CALL histclo( nid_V )
658         CALL histclo( nid_W )
659      ENDIF
660#endif
661
662   END SUBROUTINE dia_wri_state
663   !!======================================================================
664END MODULE diawri