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_c1d.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/C1D_SRC – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/C1D_SRC/diawri_c1d.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: 19.3 KB
Line 
1MODULE diawri_c1d
2   !!======================================================================
3   !!                     ***  MODULE  diawri_c1d  ***
4   !! Ocean diagnostics :  write ocean output files
5   !!=====================================================================
6   !! History :   2.0  !  2004-10  (C. Ethe)  1D Configuration
7   !!             3.0  !  2008-04  (G. Madec)  adaptation to SBC
8   !!----------------------------------------------------------------------
9#if defined key_c1d
10   !!----------------------------------------------------------------------
11   !!   'key_c1d'                                          1D Configuration
12   !!---------------------------------------------------------------------- 
13   !!   dia_wri_c1d       : create the standart NetCDF output files
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE zdf_oce         ! ocean vertical physics
18   USE sbc_oce         ! surface boundary condition: ocean
19   USE sbc_ice         ! surface boundary condition: ice
20   USE sbcmod          ! surface Boundary Codition
21   USE sbcssr          ! surface boundary condition: restauring to SSS and or SST
22   USE zdftke          ! TKE vertical mixing
23   USE zdfkpp          ! KPP vertical mixing
24   USE sol_oce         ! solver variables
25   USE phycst          ! physical constants
26   USE zdfmxl          ! mixed layer
27   USE dianam          ! build name of file (routine)
28   USE diawri
29   USE zdfddm          ! vertical  physics: double diffusion
30   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
31   USE in_out_manager  ! I/O manager
32   USE ioipsl
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC dia_wri_c1d                 ! routines called by step.F90
38
39   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT, ndex(1)   ! grid_T file
40   INTEGER, DIMENSION(jpi*jpj)     ::   ndex_hT
41   INTEGER, DIMENSION(jpi*jpj*jpk) ::   ndex_T
42
43   !! * Substitutions
44#  include "zdfddm_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/C1D 3.3 , NEMO Consortium (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50
51CONTAINS
52
53   SUBROUTINE dia_wri_c1d( kt, kindic )
54      !!---------------------------------------------------------------------
55      !!                  ***  ROUTINE dia_wri_c1d  ***
56      !!                   
57      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
58      !!      NETCDF format is used by default
59      !!
60      !! ** Method  :   At the beginning of the first time step (nit000),
61      !!      define all the NETCDF files and fields
62      !!      At each time step call histdef to compute the mean if ncessary
63      !!      Each nwrite time step, output the instantaneous or mean fields
64      !!      IF kindic <0, output of fields before the model interruption.
65      !!      IF kindic =0, time step loop
66      !!      IF kindic >0, output of fields before the time step loop
67      !!----------------------------------------------------------------------
68      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
69      INTEGER, INTENT( in ) ::   kindic  !
70      !!
71      LOGICAL ::   ll_print = .FALSE.                ! =T print and flush numout
72      CHARACTER (len=40) ::   clhstnam, clop, clmx   ! temporary names
73      INTEGER ::   inum = 11                         ! temporary logical unit
74      INTEGER ::   ji, jj, ik                        ! dummy loop indices
75      INTEGER ::   iimi, iima, ipk, it, ijmi, ijma   ! temporary integers
76      INTEGER ::   itmod                             !
77      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt   ! temporary scalars
78      REAL(wp), DIMENSION(jpi,jpj) ::   zw2d         ! temporary workspace
79      !!----------------------------------------------------------------------
80     
81      ! 0. Initialisation
82      ! -----------------
83     
84      ! local variable for debugging
85      ll_print = .FALSE.
86      ll_print = ll_print .AND. lwp
87
88      ! Define frequency of output and means
89      zdt = rdt
90      IF( nacc == 1 ) zdt = rdtmin
91      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
92      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
93      ENDIF
94#if defined key_diainstant
95      zsto = nwrite * zdt
96      clop = "inst("//TRIM(clop)//")"
97#else
98      zsto=zdt
99      clop = "ave("//TRIM(clop)//")"
100#endif
101      zout = nwrite * zdt
102      zmax = ( nitend - nit000 + 1 ) * zdt
103
104      ! Define indices of the horizontal output zoom and vertical limit storage
105      iimi = 1      ;      iima = jpi
106      ijmi = 1      ;      ijma = jpj
107      ipk = jpk
108
109      ! define time axis
110      it = kt
111      itmod = kt - nit000 + 1
112
113
114      ! 1. Define NETCDF files and fields at beginning of first time step
115      ! -----------------------------------------------------------------
116
117      IF(ll_print) WRITE(numout,*) 'dia_wri_c1d kt = ', kt, ' kindic ', kindic
118
119      IF( kt == nit000 ) THEN
120
121         ! Define the NETCDF files (one per grid)
122         
123         ! Compute julian date from starting date of the run
124         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
125         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
126         IF(lwp)WRITE(numout,*)
127         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
128            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
129         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
130                                 ' limit storage in depth = ', ipk
131
132         ! WRITE root name in date.file for use by postpro
133         CALL dia_nam( clhstnam, nwrite,' ' )
134         CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp )
135         WRITE(inum,*) clhstnam
136         CLOSE(inum)
137         
138         ! Define the T grid FILE ( nid_T )
139         
140         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
141         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
142         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
143            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
144            &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
145         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
146            &           "m", ipk, gdept_0, nz_T, "down" )
147         !                                                            ! Index of ocean points
148         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
149         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
150
151
152         ! Declare all the output fields as NETCDF variables
153
154         !                                                                                      !!! nid_T : 3D
155         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
156            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
157         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
158            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
159         !                                                                                      !!! nid_T : 2D
160         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
161            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
162         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
163            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
164         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! emp
165            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
166!!       CALL histdef( nid_T, "sorunoff", "Runoffs"                            , "Kg/m2/s",   &  ! runoffs
167!!          &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
168         CALL histdef( nid_T, "sowaflcd", "concentration/dilution water flux"  , "kg/m2/s",   &  ! emps
169            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
170         CALL histdef( nid_T, "sosalflx", "Surface Salt Flux"                  , "Kg/m2/s",   &  ! emps * sn
171            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
172         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qsr + qns
173            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
174         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
175            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
176         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
177            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
178#if defined key_zdfkpp
179         CALL histdef( nid_T, "sokppekd", "Ekman depth                     "   , "m"      ,   &  ! sokppekd
180            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
181         CALL histdef( nid_T, "sokppbld", "Boundary Layer Depth            "   , "m"      ,   &  ! sokppbld
182            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
183#endif
184         CALL histdef( nid_T, "somxlavt", "AVT : bottom of the mixed layer    ", "m"      ,   &  ! avt_mxl
185            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
186         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
187            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
188         CALL histdef( nid_T, "soicecov", "Ice Fraction"                          , "[0,1]"  ,   &  ! fr_i
189            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
190         IF( ln_ssr ) THEN
191            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
192               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
193            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
194               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
195            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
196               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
197         ENDIF
198         clmx ="l_max(only(x))"    ! max index on a period
199         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
200            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
201         !                                                                                      !!! nid_U : 3D
202         CALL histdef( nid_T, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
203            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
204#if defined key_diaeiv
205         CALL histdef( nid_T, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
206            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
207#endif
208         !                                                                                      !!! nid_U : 2D
209         CALL histdef( nid_T, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
210            &          jpi, jpj, nh_T, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
211
212         !                                                                                      !!! nid_V : 3D
213         CALL histdef( nid_T, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
214            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
215#if defined key_diaeiv
216         CALL histdef( nid_T, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
217            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
218#endif
219         !                                                                                      !!! nid_V : 2D
220         CALL histdef( nid_T, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
221            &          jpi, jpj, nh_T, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
222#if defined key_zdfkpp
223         CALL histdef( nid_T, "vokpprig", " Gradient Richardson Number"        ,  "-"      ,   &  ! rig
224            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
225         CALL histdef( nid_T, "vokpprib", " Bulk Richardson Number    "        ,  "-"      ,   &   ! rib
226            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
227         CALL histdef( nid_T, "vokppbsf", " Buoyancy forcing          "        , "N/m2"    ,   &  ! sokppbsf
228            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
229         CALL histdef( nid_T, "vokppmol", "Moning Obukhov length scale     "   , "m"       ,   &  ! sokppmol
230            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
231#endif
232         !
233         CALL histdef( nid_T, "voeosbn2", "Brunt-Vaisala Frequency"             , "m2/s2"  ,   &  ! rn2
234            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
235         CALL histdef( nid_T, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
236            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
237         CALL histdef( nid_T, "votkeavm", "Vertical Eddy Viscosity",             "m2/s"   ,   &  ! avmu
238            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
239         !
240         IF( lk_zdfddm ) THEN
241            CALL histdef( nid_T,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
242               &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
243         ENDIF
244
245         CALL histend( nid_T, snc4set )
246
247         IF(lwp) WRITE(numout,*)
248         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
249         IF(ll_print) CALL FLUSH(numout )
250
251      ENDIF
252
253      ! 2. Start writing data
254      ! ---------------------
255
256      ! ndex(1) est utilise ssi l'avant dernier argument est diffferent de
257      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
258      ! donne le nombre d'elements, et ndex la liste des indices a sortir
259
260      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
261         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
262         WRITE(numout,*) '~~~~~~ '
263      ENDIF
264
265      ! Write fields on T grid
266      CALL histwrite( nid_T, "votemper", it, tn            , ndim_T , ndex_T  )   ! temperature
267      CALL histwrite( nid_T, "vosaline", it, sn            , ndim_T , ndex_T  )   ! salinity
268      CALL histwrite( nid_T, "sosstsst", it, tn(:,:,1)     , ndim_hT, ndex_hT )   ! sea surface temperature
269      CALL histwrite( nid_T, "sosaline", it, sn(:,:,1)     , ndim_hT, ndex_hT )   ! sea surface salinity
270      CALL histwrite( nid_T, "sowaflup", it, emp           , ndim_hT, ndex_hT )   ! upward water flux
271!!    CALL histwrite( nid_T, "sorunoff", it, runoff        , ndim_hT, ndex_hT )   ! runoff
272      CALL histwrite( nid_T, "sowaflcd", it, emps          , ndim_hT, ndex_hT )   ! c/d water flux
273      zw2d(:,:) = emps(:,:) * sn(:,:,1) * tmask(:,:,1)
274      CALL histwrite( nid_T, "sosalflx", it, zw2d          , ndim_hT, ndex_hT )   ! c/d salt flux
275      CALL histwrite( nid_T, "sohefldo", it, qsr + qns     , ndim_hT, ndex_hT )   ! total heat flux
276      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
277      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
278#if defined key_zdfkpp
279      CALL histwrite( nid_T, "sokppekd", it, ekdp          , ndim_hT, ndex_hT )   ! Ekman depht
280      CALL histwrite( nid_T, "sokppbld", it, hkpp          , ndim_hT, ndex_hT )   ! boundary layer depth
281#endif 
282      ! store the vertical eddy diffusivity coef. at the bottom of the mixed layer
283      DO jj = 1, jpj
284         DO ji = 1, jpi
285            ik = nmln(ji,jj)
286            zw2d(ji,jj) = avt(ji,jj,ik) * tmask(ji,jj,1)
287         END DO
288      END DO
289      CALL histwrite( nid_T, "somxlavt", it, zw2d          , ndim_hT, ndex_hT )   ! Kz at bottom of mixed layer
290      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
291      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction
292      IF( ln_ssr ) THEN
293         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
294         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
295         zw2d(:,:) = erp(:,:) * sn(:,:,1) * tmask(:,:,1)
296         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
297      ENDIF
298      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
299      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
300      CALL histwrite( nid_T, "vozocrtx", it, un            , ndim_T , ndex_T )    ! i-current
301      CALL histwrite( nid_T, "sozotaux", it, utau          , ndim_hT, ndex_hT )   ! i-wind stress
302      CALL histwrite( nid_T, "vomecrty", it, vn            , ndim_T , ndex_T  )   ! j-current
303      CALL histwrite( nid_T, "sometauy", it, vtau          , ndim_hT, ndex_hT )   ! j-wind stress
304#if defined key_zdfkpp
305      CALL histwrite( nid_T, "vokpprig", it, rig           , ndim_T , ndex_T )    ! gradient Richardson number
306      CALL histwrite( nid_T, "vokpprib", it, rib           , ndim_T , ndex_T )    ! bulk Richardson number
307      CALL histwrite( nid_T, "vokppbsf", it, buof          , ndim_T , ndex_T )    ! buoyancy forcing
308      CALL histwrite( nid_T, "vokppmol", it, mols          , ndim_T , ndex_T )    ! Moning-Obukov length scale
309#endif
310      CALL histwrite( nid_T, "voeosbn2", it, rn2           , ndim_T , ndex_T )    ! Brunt-Vaisala Frequency
311      CALL histwrite( nid_T, "votkeavt", it, avt           , ndim_T , ndex_T )    ! T vert. eddy diff. coef.
312      CALL histwrite( nid_T, "votkeavm", it, avmu          , ndim_T , ndex_T )    ! T vert. eddy visc. coef.
313      IF( lk_zdfddm ) THEN
314         CALL histwrite( nid_T, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
315      ENDIF
316
317      !  Create an output files (output.abort.nc) if S < 0 or u > 20 m/s
318      IF( kindic < 0 )   CALL dia_wri_state( 'output.abort', kt )
319
320      ! 3. Close all files
321      ! ---------------------------------------
322      IF( kt == nitend .OR. kindic < 0 )   CALL histclo( nid_T )
323      !
324   END SUBROUTINE dia_wri_c1d
325
326#else
327   !!----------------------------------------------------------------------
328   !!   Default key                                     NO 1D Configuration
329   !!----------------------------------------------------------------------
330CONTAINS
331   SUBROUTINE dia_wri_c1d ( kt, kindic )      ! dummy routine
332      WRITE(*,*) 'dia_wri_c1d: You should not have seen this print! error?', kt, kindic
333   END SUBROUTINE dia_wri_c1d
334#endif
335
336   !!======================================================================
337END MODULE diawri_c1d
Note: See TracBrowser for help on using the repository browser.