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

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

source: trunk/NEMO/OPA_SRC/DIA/diawri.F90 @ 1756

Last change on this file since 1756 was 1756, checked in by smasson, 14 years ago

implement AR5 diagnostics, see ticket:610

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