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

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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 2364

Last change on this file since 2364 was 2364, checked in by acc, 13 years ago

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