source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/C1D_SRC/diawri_c1d.F90 @ 2382

Last change on this file since 2382 was 2382, checked in by gm, 10 years ago

v3.3beta: C1D - bug correction to compile with key_c1d

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