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 NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diawri.F90 @ 11822

Last change on this file since 11822 was 11822, checked in by acc, 4 years ago

Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Sette tested updates to branch to align with trunk changes between 10721 and 11740. Sette tests are passing but results differ from branch before these changes (except for GYRE_PISCES and VORTEX) and branch results already differed from trunk because of algorithmic fixes. Will need more checks to confirm correctness.

  • Property svn:keywords set to Id
File size: 49.7 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   !!            3.7  ! 2014-01  (G. Madec) remove eddy induced velocity from no-IOM output
20   !!                 !                     change name of output variables in dia_wri_state
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   dia_wri       : create the standart output files
25   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
26   !!----------------------------------------------------------------------
27   USE oce            ! ocean dynamics and tracers
28   USE dom_oce        ! ocean space and time domain
29   USE phycst         ! physical constants
30   USE dianam         ! build name of file (routine)
31   USE diahth         ! thermocline diagnostics
32   USE dynadv   , ONLY: ln_dynadv_vec
33   USE icb_oce        ! Icebergs
34   USE icbdia         ! Iceberg budgets
35   USE ldftra         ! lateral physics: eddy diffusivity coef.
36   USE ldfdyn         ! lateral physics: eddy viscosity   coef.
37   USE sbc_oce        ! Surface boundary condition: ocean fields
38   USE sbc_ice        ! Surface boundary condition: ice fields
39   USE sbcssr         ! restoring term toward SST/SSS climatology
40   USE sbcwave        ! wave parameters
41   USE wet_dry        ! wetting and drying
42   USE zdf_oce        ! ocean vertical physics
43   USE zdfdrg         ! ocean vertical physics: top/bottom friction
44   USE zdfmxl         ! mixed layer
45   !
46   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
47   USE in_out_manager ! I/O manager
48   USE diatmb         ! Top,middle,bottom output
49   USE dia25h         ! 25h Mean output
50   USE iom            !
51   USE ioipsl         !
52
53#if defined key_si3
54   USE ice 
55   USE icewri 
56#endif
57   USE lib_mpp         ! MPP library
58   USE timing          ! preformance summary
59   USE diu_bulk        ! diurnal warm layer
60   USE diu_coolskin    ! Cool skin
61
62   IMPLICIT NONE
63   PRIVATE
64
65   PUBLIC   dia_wri                 ! routines called by step.F90
66   PUBLIC   dia_wri_state
67   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
68
69   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
70   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
71   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
72   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
73   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
74   INTEGER ::   ndex(1)                              ! ???
75   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
76   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
77   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
78
79   !! * Substitutions
80#  include "vectopt_loop_substitute.h90"
81   !!----------------------------------------------------------------------
82   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
83   !! $Id$
84   !! Software governed by the CeCILL license (see ./LICENSE)
85   !!----------------------------------------------------------------------
86CONTAINS
87
88#if defined key_iomput
89   !!----------------------------------------------------------------------
90   !!   'key_iomput'                                        use IOM library
91   !!----------------------------------------------------------------------
92   INTEGER FUNCTION dia_wri_alloc()
93      !
94      dia_wri_alloc = 0
95      !
96   END FUNCTION dia_wri_alloc
97
98   
99   SUBROUTINE dia_wri( kt, Kmm )
100      !!---------------------------------------------------------------------
101      !!                  ***  ROUTINE dia_wri  ***
102      !!                   
103      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
104      !!      NETCDF format is used by default
105      !!
106      !! ** Method  :  use iom_put
107      !!----------------------------------------------------------------------
108      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
109      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index
110      !!
111      INTEGER ::   ji, jj, jk       ! dummy loop indices
112      INTEGER ::   ikbot            ! local integer
113      REAL(wp)::   zztmp , zztmpx   ! local scalar
114      REAL(wp)::   zztmp2, zztmpy   !   -      -
115      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
116      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
117      !!----------------------------------------------------------------------
118      !
119      IF( ln_timing )   CALL timing_start('dia_wri')
120      !
121      ! Output the initial state and forcings
122      IF( ninist == 1 ) THEN                       
123         CALL dia_wri_state( Kmm, 'output.init' )
124         ninist = 0
125      ENDIF
126
127      ! Output of initial vertical scale factor
128      CALL iom_put("e3t_0", e3t_0(:,:,:) )
129      CALL iom_put("e3u_0", e3u_0(:,:,:) )
130      CALL iom_put("e3v_0", e3v_0(:,:,:) )
131      !
132      CALL iom_put( "e3t" , e3t(:,:,:,Kmm) )
133      CALL iom_put( "e3u" , e3u(:,:,:,Kmm) )
134      CALL iom_put( "e3v" , e3v(:,:,:,Kmm) )
135      CALL iom_put( "e3w" , e3w(:,:,:,Kmm) )
136      IF( iom_use("e3tdef") )   &
137         CALL iom_put( "e3tdef"  , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
138
139      IF( ll_wd ) THEN
140         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
141      ELSE
142         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height
143      ENDIF
144
145      IF( iom_use("wetdep") )   &                  ! wet depth
146         CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) )
147     
148      CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) )    ! 3D temperature
149      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature
150      IF ( iom_use("sbt") ) THEN
151         DO jj = 1, jpj
152            DO ji = 1, jpi
153               ikbot = mbkt(ji,jj)
154               z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
155            END DO
156         END DO
157         CALL iom_put( "sbt", z2d )                ! bottom temperature
158      ENDIF
159     
160      CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) )    ! 3D salinity
161      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity
162      IF ( iom_use("sbs") ) THEN
163         DO jj = 1, jpj
164            DO ji = 1, jpi
165               ikbot = mbkt(ji,jj)
166               z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
167            END DO
168         END DO
169         CALL iom_put( "sbs", z2d )                ! bottom salinity
170      ENDIF
171
172      IF ( iom_use("taubot") ) THEN                ! bottom stress
173         zztmp = rau0 * 0.25
174         z2d(:,:) = 0._wp
175         DO jj = 2, jpjm1
176            DO ji = fs_2, fs_jpim1   ! vector opt.
177               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   &
178                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   &
179                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   &
180                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2
181               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
182               !
183            END DO
184         END DO
185         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
186         CALL iom_put( "taubot", z2d )           
187      ENDIF
188         
189      CALL iom_put( "uoce", uu(:,:,:,Kmm) )            ! 3D i-current
190      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current
191      IF ( iom_use("sbu") ) THEN
192         DO jj = 1, jpj
193            DO ji = 1, jpi
194               ikbot = mbku(ji,jj)
195               z2d(ji,jj) = uu(ji,jj,ikbot,Kmm)
196            END DO
197         END DO
198         CALL iom_put( "sbu", z2d )                ! bottom i-current
199      ENDIF
200     
201      CALL iom_put( "voce", vv(:,:,:,Kmm) )            ! 3D j-current
202      CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current
203      IF ( iom_use("sbv") ) THEN
204         DO jj = 1, jpj
205            DO ji = 1, jpi
206               ikbot = mbkv(ji,jj)
207               z2d(ji,jj) = vv(ji,jj,ikbot,Kmm)
208            END DO
209         END DO
210         CALL iom_put( "sbv", z2d )                ! bottom j-current
211      ENDIF
212
213      IF( ln_zad_Aimp ) ww = ww + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
214      !
215      CALL iom_put( "woce", ww )                   ! vertical velocity
216      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
217         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
218         z2d(:,:) = rau0 * e1e2t(:,:)
219         DO jk = 1, jpk
220            z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:)
221         END DO
222         CALL iom_put( "w_masstr" , z3d ) 
223         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
224      ENDIF
225      !
226      IF( ln_zad_Aimp ) ww = ww - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output
227
228      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
229      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef.
230      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef.
231
232      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
233      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
234
235      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
236         DO jj = 2, jpjm1                                    ! sst gradient
237            DO ji = fs_2, fs_jpim1   ! vector opt.
238               zztmp  = ts(ji,jj,1,jp_tem,Kmm)
239               zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
240               zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1)
241               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
242                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
243            END DO
244         END DO
245         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
246         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient
247         z2d(:,:) = SQRT( z2d(:,:) )
248         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient
249      ENDIF
250         
251      ! heat and salt contents
252      IF( iom_use("heatc") ) THEN
253         z2d(:,:)  = 0._wp 
254         DO jk = 1, jpkm1
255            DO jj = 1, jpj
256               DO ji = 1, jpi
257                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
258               END DO
259            END DO
260         END DO
261         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
262      ENDIF
263
264      IF( iom_use("saltc") ) THEN
265         z2d(:,:)  = 0._wp 
266         DO jk = 1, jpkm1
267            DO jj = 1, jpj
268               DO ji = 1, jpi
269                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
270               END DO
271            END DO
272         END DO
273         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
274      ENDIF
275      !
276      IF ( iom_use("eken") ) THEN
277         z3d(:,:,jpk) = 0._wp 
278         DO jk = 1, jpkm1
279            DO jj = 2, jpjm1
280               DO ji = fs_2, fs_jpim1   ! vector opt.
281                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
282                  z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   &
283                     &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   &
284                     &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   &
285                     &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   )
286               END DO
287            END DO
288         END DO
289         CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
290         CALL iom_put( "eken", z3d )                 ! kinetic energy
291      ENDIF
292      !
293      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence
294      !
295      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
296         z3d(:,:,jpk) = 0.e0
297         z2d(:,:) = 0.e0
298         DO jk = 1, jpkm1
299            z3d(:,:,jk) = rau0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
300            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
301         END DO
302         CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction
303         CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum
304      ENDIF
305     
306      IF( iom_use("u_heattr") ) THEN
307         z2d(:,:) = 0._wp 
308         DO jk = 1, jpkm1
309            DO jj = 2, jpjm1
310               DO ji = fs_2, fs_jpim1   ! vector opt.
311                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
312               END DO
313            END DO
314         END DO
315         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
316         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction
317      ENDIF
318
319      IF( iom_use("u_salttr") ) THEN
320         z2d(:,:) = 0.e0 
321         DO jk = 1, jpkm1
322            DO jj = 2, jpjm1
323               DO ji = fs_2, fs_jpim1   ! vector opt.
324                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
325               END DO
326            END DO
327         END DO
328         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
329         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction
330      ENDIF
331
332     
333      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
334         z3d(:,:,jpk) = 0.e0
335         DO jk = 1, jpkm1
336            z3d(:,:,jk) = rau0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
337         END DO
338         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
339      ENDIF
340     
341      IF( iom_use("v_heattr") ) THEN
342         z2d(:,:) = 0.e0 
343         DO jk = 1, jpkm1
344            DO jj = 2, jpjm1
345               DO ji = fs_2, fs_jpim1   ! vector opt.
346                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
347               END DO
348            END DO
349         END DO
350         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
351         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction
352      ENDIF
353
354      IF( iom_use("v_salttr") ) THEN
355         z2d(:,:) = 0._wp 
356         DO jk = 1, jpkm1
357            DO jj = 2, jpjm1
358               DO ji = fs_2, fs_jpim1   ! vector opt.
359                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
360               END DO
361            END DO
362         END DO
363         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
364         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction
365      ENDIF
366
367      IF( iom_use("tosmint") ) THEN
368         z2d(:,:) = 0._wp
369         DO jk = 1, jpkm1
370            DO jj = 2, jpjm1
371               DO ji = fs_2, fs_jpim1   ! vector opt.
372                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm)
373               END DO
374            END DO
375         END DO
376         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
377         CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature
378      ENDIF
379      IF( iom_use("somint") ) THEN
380         z2d(:,:)=0._wp
381         DO jk = 1, jpkm1
382            DO jj = 2, jpjm1
383               DO ji = fs_2, fs_jpim1   ! vector opt.
384                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
385               END DO
386            END DO
387         END DO
388         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
389         CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity
390      ENDIF
391
392      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2)
393      !
394
395      IF (ln_diatmb)   CALL dia_tmb( Kmm )            ! tmb values
396         
397      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging
398
399      IF( ln_timing )   CALL timing_stop('dia_wri')
400      !
401   END SUBROUTINE dia_wri
402
403#else
404   !!----------------------------------------------------------------------
405   !!   Default option                                  use IOIPSL  library
406   !!----------------------------------------------------------------------
407
408   INTEGER FUNCTION dia_wri_alloc()
409      !!----------------------------------------------------------------------
410      INTEGER, DIMENSION(2) :: ierr
411      !!----------------------------------------------------------------------
412      ierr = 0
413      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
414         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
415         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
416         !
417      dia_wri_alloc = MAXVAL(ierr)
418      CALL mpp_sum( 'diawri', dia_wri_alloc )
419      !
420   END FUNCTION dia_wri_alloc
421
422   
423   SUBROUTINE dia_wri( kt, Kmm )
424      !!---------------------------------------------------------------------
425      !!                  ***  ROUTINE dia_wri  ***
426      !!                   
427      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
428      !!      NETCDF format is used by default
429      !!
430      !! ** Method  :   At the beginning of the first time step (nit000),
431      !!      define all the NETCDF files and fields
432      !!      At each time step call histdef to compute the mean if ncessary
433      !!      Each nn_write time step, output the instantaneous or mean fields
434      !!----------------------------------------------------------------------
435      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
436      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index
437      !
438      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
439      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
440      INTEGER  ::   inum = 11                                ! temporary logical unit
441      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
442      INTEGER  ::   ierr                                     ! error code return from allocation
443      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
444      INTEGER  ::   jn, ierror                               ! local integers
445      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
446      !
447      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
448      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
449      !!----------------------------------------------------------------------
450      !
451      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
452         CALL dia_wri_state( Kmm, 'output.init' )
453         ninist = 0
454      ENDIF
455      !
456      IF( nn_write == -1 )   RETURN   ! we will never do any output
457      !
458      IF( ln_timing )   CALL timing_start('dia_wri')
459      !
460      ! 0. Initialisation
461      ! -----------------
462
463      ll_print = .FALSE.                  ! local variable for debugging
464      ll_print = ll_print .AND. lwp
465
466      ! Define frequency of output and means
467      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
468#if defined key_diainstant
469      zsto = nn_write * rdt
470      clop = "inst("//TRIM(clop)//")"
471#else
472      zsto=rdt
473      clop = "ave("//TRIM(clop)//")"
474#endif
475      zout = nn_write * rdt
476      zmax = ( nitend - nit000 + 1 ) * rdt
477
478      ! Define indices of the horizontal output zoom and vertical limit storage
479      iimi = 1      ;      iima = jpi
480      ijmi = 1      ;      ijma = jpj
481      ipk = jpk
482
483      ! define time axis
484      it = kt
485      itmod = kt - nit000 + 1
486
487
488      ! 1. Define NETCDF files and fields at beginning of first time step
489      ! -----------------------------------------------------------------
490
491      IF( kt == nit000 ) THEN
492
493         ! Define the NETCDF files (one per grid)
494
495         ! Compute julian date from starting date of the run
496         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
497         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
498         IF(lwp)WRITE(numout,*)
499         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
500            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
501         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
502                                 ' limit storage in depth = ', ipk
503
504         ! WRITE root name in date.file for use by postpro
505         IF(lwp) THEN
506            CALL dia_nam( clhstnam, nn_write,' ' )
507            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
508            WRITE(inum,*) clhstnam
509            CLOSE(inum)
510         ENDIF
511
512         ! Define the T grid FILE ( nid_T )
513
514         CALL dia_nam( clhstnam, nn_write, 'grid_T' )
515         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
516         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
517            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
518            &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
519         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
520            &           "m", ipk, gdept_1d, nz_T, "down" )
521         !                                                            ! Index of ocean points
522         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
523         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
524         !
525         IF( ln_icebergs ) THEN
526            !
527            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
528            !! that routine is called from nemogcm, so do it here immediately before its needed
529            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
530            CALL mpp_sum( 'diawri', ierror )
531            IF( ierror /= 0 ) THEN
532               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
533               RETURN
534            ENDIF
535            !
536            !! iceberg vertical coordinate is class number
537            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
538               &           "number", nclasses, class_num, nb_T )
539            !
540            !! each class just needs the surface index pattern
541            ndim_bT = 3
542            DO jn = 1,nclasses
543               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
544            ENDDO
545            !
546         ENDIF
547
548         ! Define the U grid FILE ( nid_U )
549
550         CALL dia_nam( clhstnam, nn_write, 'grid_U' )
551         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
552         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
553            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
554            &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
555         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
556            &           "m", ipk, gdept_1d, nz_U, "down" )
557         !                                                            ! Index of ocean points
558         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
559         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
560
561         ! Define the V grid FILE ( nid_V )
562
563         CALL dia_nam( clhstnam, nn_write, 'grid_V' )                   ! filename
564         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
565         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
566            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
567            &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
568         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
569            &          "m", ipk, gdept_1d, nz_V, "down" )
570         !                                                            ! Index of ocean points
571         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
572         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
573
574         ! Define the W grid FILE ( nid_W )
575
576         CALL dia_nam( clhstnam, nn_write, 'grid_W' )                   ! filename
577         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
578         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
579            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
580            &          nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
581         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
582            &          "m", ipk, gdepw_1d, nz_W, "down" )
583
584
585         ! Declare all the output fields as NETCDF variables
586
587         !                                                                                      !!! nid_T : 3D
588         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
589            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
590         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
591            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
592         IF(  .NOT.ln_linssh  ) THEN
593            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t(:,:,:,Kmm)
594            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
595            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t(:,:,:,Kmm)
596            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
597            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t(:,:,:,Kmm)
598            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
599         ENDIF
600         !                                                                                      !!! nid_T : 2D
601         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
602            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
603         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
604            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
605         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
606            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
607         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
608            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
609         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
610            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
611         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
612            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
613         IF(  ln_linssh  ) THEN
614            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * ts(:,:,1,jp_tem,Kmm)
615            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
616            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
617            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * ts(:,:,1,jp_sal,Kmm)
618            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
619            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
620         ENDIF
621         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
622            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
623         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
624            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
625         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
626            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
627         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
628            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
629         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
630            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
631         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
632            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
633!
634         IF( ln_icebergs ) THEN
635            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
636               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
637            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
638               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
639            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
640               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
641            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
642               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
643            IF( ln_bergdia ) THEN
644               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
645                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
646               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
647                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
648               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
649                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
650               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
651                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
652               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
653                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
654               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
655                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
656               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
657                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
658               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
659                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
660               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
661                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
662               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
663                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
664            ENDIF
665         ENDIF
666
667         IF( ln_ssr ) THEN
668            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
669               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
670            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
671               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
672            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
673               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
674         ENDIF
675       
676         clmx ="l_max(only(x))"    ! max index on a period
677!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
678!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
679#if defined key_diahth
680         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
681            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
682         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
683            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
684         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
685            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
686         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "J/m2"   ,   & ! htc3
687            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
688#endif
689
690         CALL histend( nid_T, snc4chunks=snc4set )
691
692         !                                                                                      !!! nid_U : 3D
693         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! uu(:,:,:,Kmm)
694            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
695         IF( ln_wave .AND. ln_sdw) THEN
696            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd
697               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
698         ENDIF
699         !                                                                                      !!! nid_U : 2D
700         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
701            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
702
703         CALL histend( nid_U, snc4chunks=snc4set )
704
705         !                                                                                      !!! nid_V : 3D
706         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vv(:,:,:,Kmm)
707            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
708         IF( ln_wave .AND. ln_sdw) THEN
709            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd
710               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
711         ENDIF
712         !                                                                                      !!! nid_V : 2D
713         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
714            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
715
716         CALL histend( nid_V, snc4chunks=snc4set )
717
718         !                                                                                      !!! nid_W : 3D
719         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! ww
720            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
721         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
722            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
723         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avm
724            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
725
726         IF( ln_zdfddm ) THEN
727            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
728               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
729         ENDIF
730         
731         IF( ln_wave .AND. ln_sdw) THEN
732            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd
733               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
734         ENDIF
735         !                                                                                      !!! nid_W : 2D
736         CALL histend( nid_W, snc4chunks=snc4set )
737
738         IF(lwp) WRITE(numout,*)
739         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
740         IF(ll_print) CALL FLUSH(numout )
741
742      ENDIF
743
744      ! 2. Start writing data
745      ! ---------------------
746
747      ! ndex(1) est utilise ssi l'avant dernier argument est different de
748      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
749      ! donne le nombre d'elements, et ndex la liste des indices a sortir
750
751      IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
752         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
753         WRITE(numout,*) '~~~~~~ '
754      ENDIF
755
756      IF( .NOT.ln_linssh ) THEN
757         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! heat content
758         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! salt content
759         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface heat content
760         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity content
761      ELSE
762         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature
763         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T  )   ! salinity
764         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT )   ! sea surface temperature
765         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity
766      ENDIF
767      IF( .NOT.ln_linssh ) THEN
768         zw3d(:,:,:) = ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
769         CALL histwrite( nid_T, "vovvle3t", it, e3t (:,:,:,Kmm) , ndim_T , ndex_T  )   ! level thickness
770         CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T  )   ! t-point depth
771         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
772      ENDIF
773      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)          , ndim_hT, ndex_hT )   ! sea surface height
774      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
775      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
776      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
777                                                                                  ! (includes virtual salt flux beneath ice
778                                                                                  ! in linear free surface case)
779      IF( ln_linssh ) THEN
780         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm)
781         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
782         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm)
783         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
784      ENDIF
785      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
786      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
787      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
788      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
789      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
790      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
791!
792      IF( ln_icebergs ) THEN
793         !
794         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
795         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
796         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
797         !
798         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
799         !
800         IF( ln_bergdia ) THEN
801            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
802            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
803            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
804            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
805            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
806            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
807            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
808            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
809            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
810            !
811            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
812         ENDIF
813      ENDIF
814
815      IF( ln_ssr ) THEN
816         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
817         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
818         zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
819         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
820      ENDIF
821!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
822!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
823
824#if defined key_diahth
825      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
826      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
827      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
828      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
829#endif
830
831      CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm)            , ndim_U , ndex_U )    ! i-current
832      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
833
834      CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm)            , ndim_V , ndex_V  )   ! j-current
835      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
836
837      IF( ln_zad_Aimp ) THEN
838         CALL histwrite( nid_W, "vovecrtz", it, ww + wi     , ndim_T, ndex_T )    ! vert. current
839      ELSE
840         CALL histwrite( nid_W, "vovecrtz", it, ww          , ndim_T, ndex_T )    ! vert. current
841      ENDIF
842      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
843      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
844      IF( ln_zdfddm ) THEN
845         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef.
846      ENDIF
847
848      IF( ln_wave .AND. ln_sdw ) THEN
849         CALL histwrite( nid_U, "sdzocrtx", it, usd         , ndim_U , ndex_U )    ! i-StokesDrift-current
850         CALL histwrite( nid_V, "sdmecrty", it, vsd         , ndim_V , ndex_V )    ! j-StokesDrift-current
851         CALL histwrite( nid_W, "sdvecrtz", it, wsd         , ndim_T , ndex_T )    ! StokesDrift vert. current
852      ENDIF
853
854      ! 3. Close all files
855      ! ---------------------------------------
856      IF( kt == nitend ) THEN
857         CALL histclo( nid_T )
858         CALL histclo( nid_U )
859         CALL histclo( nid_V )
860         CALL histclo( nid_W )
861      ENDIF
862      !
863      IF( ln_timing )   CALL timing_stop('dia_wri')
864      !
865   END SUBROUTINE dia_wri
866#endif
867
868   SUBROUTINE dia_wri_state( Kmm, cdfile_name )
869      !!---------------------------------------------------------------------
870      !!                 ***  ROUTINE dia_wri_state  ***
871      !!       
872      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
873      !!      the instantaneous ocean state and forcing fields.
874      !!        Used to find errors in the initial state or save the last
875      !!      ocean state in case of abnormal end of a simulation
876      !!
877      !! ** Method  :   NetCDF files using ioipsl
878      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
879      !!      File 'output.abort.nc' is created in case of abnormal job end
880      !!----------------------------------------------------------------------
881      INTEGER           , INTENT( in ) ::   Kmm              ! time level index
882      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
883      !!
884      INTEGER :: inum
885      !!----------------------------------------------------------------------
886      !
887      IF(lwp) WRITE(numout,*)
888      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
889      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
890      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
891
892#if defined key_si3
893     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
894#else
895     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
896#endif
897
898      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature
899      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity
900      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height
901      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity
902      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity
903      IF( ln_zad_Aimp ) THEN
904         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi        )    ! now k-velocity
905      ELSE
906         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity
907      ENDIF
908      IF( ALLOCATED(ahtu) ) THEN
909         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
910         CALL iom_rstput( 0, 0, inum,  'ahtv', ahtv              )    ! aht at v-point
911      ENDIF
912      IF( ALLOCATED(ahmt) ) THEN
913         CALL iom_rstput( 0, 0, inum,  'ahmt', ahmt              )    ! ahmt at u-point
914         CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
915      ENDIF
916      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
917      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
918      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
919      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
920      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
921      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
922      IF(  .NOT.ln_linssh  ) THEN             
923         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth
924         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm)          )    !  T-cell thickness 
925      END IF
926      IF( ln_wave .AND. ln_sdw ) THEN
927         CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
928         CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd            )    ! now StokesDrift j-velocity
929         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
930      ENDIF
931 
932#if defined key_si3
933      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
934         CALL ice_wri_state( inum )
935      ENDIF
936#endif
937      !
938      CALL iom_close( inum )
939      !
940   END SUBROUTINE dia_wri_state
941
942   !!======================================================================
943END MODULE diawri
Note: See TracBrowser for help on using the repository browser.