source: trunk/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90 @ 3564

Last change on this file since 3564 was 3564, checked in by rblod, 8 years ago

fix output average with EVP, see ticket #908

  • Property svn:keywords set to Id
File size: 17.5 KB
Line 
1MODULE limwri_2
2   !!======================================================================
3   !!                     ***  MODULE  limwri_2  ***
4   !!         Ice diagnostics :  write ice output files
5   !!======================================================================
6   !! history :  2.0  ! 2003-08  (C. Ethe)      original code
7   !!            2.0  ! 2004-10  (C. Ethe )     1D configuration
8   !!             -   ! 2009-06  (B. Lemaire )  iom_put + lim_wri_state_2
9   !!-------------------------------------------------------------------
10#if defined key_lim2
11   !!----------------------------------------------------------------------
12   !!   'key_lim2'                                    LIM 2.0 sea-ice model
13   !!----------------------------------------------------------------------
14   !!----------------------------------------------------------------------
15   !!   lim_wri_2      : write of the diagnostics variables in ouput file
16   !!   lim_wri_init_2 : initialization and namelist read
17   !!   lim_wri_state_2 : write for initial state or/and abandon:
18   !!                     > output.init.nc (if ninist = 1 in namelist)
19   !!                     > output.abort.nc
20   !!----------------------------------------------------------------------
21   USE phycst
22   USE dom_oce
23   USE sbc_oce
24   USE sbc_ice
25   USE dom_ice_2
26   USE ice_2
27
28   USE dianam          ! build name of file (routine)
29   USE lbclnk
30   USE in_out_manager
31   USE lib_mpp         ! MPP library
32   USE wrk_nemo        ! work arrays
33   USE iom
34   USE ioipsl
35
36   IMPLICIT NONE
37   PRIVATE
38
39#if ! defined key_iomput
40   PUBLIC   lim_wri_2         ! called by sbc_ice_lim_2
41#endif
42   PUBLIC   lim_wri_state_2   ! called by dia_wri_state
43   PUBLIC   lim_wri_alloc_2   ! called by nemogcm.F90
44
45   INTEGER, PARAMETER                       ::   jpnoumax = 40   ! maximum number of variable for ice output
46   INTEGER                                  ::   noumef          ! number of fields
47   REAL(wp)           , DIMENSION(jpnoumax) ::   cmulti ,     &  ! multiplicative constant
48      &                                          cadd            ! additive constant
49   CHARACTER(len = 35), DIMENSION(jpnoumax) ::   titn            ! title of the field
50   CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::   nam             ! name of the field
51   CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::   uni             ! unit of the field
52   INTEGER            , DIMENSION(jpnoumax) ::   nc              ! switch for saving field ( = 1 ) or not ( = 0 )
53
54   INTEGER ::   nice, nhorid, ndim, niter, ndepid       ! ????
55   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndex51   ! ????
56
57   REAL(wp) ::   epsi16 = 1.e-16_wp   ! constant values
58   REAL(wp) ::   zzero  = 0._wp       !     -      -
59   REAL(wp) ::   zone   = 1._wp       !     -      -
60
61   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zcmo      ! Workspace array for netcdf writer.
62
63
64   !! * Substitutions
65#   include "vectopt_loop_substitute.h90"
66   !!----------------------------------------------------------------------
67   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
68   !! $Id$
69   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
70   !!----------------------------------------------------------------------
71CONTAINS
72
73   INTEGER FUNCTION lim_wri_alloc_2()
74      !!-------------------------------------------------------------------
75      !!                  ***   ROUTINE lim_wri_alloc_2  ***
76      !!-------------------------------------------------------------------
77      ALLOCATE( ndex51(jpij), zcmo(jpi,jpj,jpnoumax), STAT=lim_wri_alloc_2)
78      !
79      IF( lk_mpp               )   CALL mpp_sum ( lim_wri_alloc_2 )
80      IF( lim_wri_alloc_2 /= 0 )   CALL ctl_warn('lim_wri_alloc_2: failed to allocate array ndex51')
81      !
82   END FUNCTION lim_wri_alloc_2
83
84
85#if ! defined key_iomput
86# if defined key_dimgout
87   !!----------------------------------------------------------------------
88   !!   'key_dimgout'                                    Direct Access file
89   !!----------------------------------------------------------------------
90# include "limwri_dimg_2.h90"
91# else
92   SUBROUTINE lim_wri_2( kt )
93      !!-------------------------------------------------------------------
94      !!                    ***   ROUTINE lim_wri_2  ***
95      !!               
96      !! ** Purpose :   write the sea-ice output file in NetCDF
97      !!
98      !! ** Method  :   computes the average of some variables and write
99      !!      it in the NetCDF ouput files
100      !!      CAUTION: the sea-ice time-step must be an integer fraction
101      !!      of a day
102      !!-------------------------------------------------------------------
103      INTEGER, INTENT(in) ::   kt     ! number of iteration
104      !!
105      INTEGER  ::   ji, jj, jf                      ! dummy loop indices
106      CHARACTER(len = 80)  ::   clhstnam, clop
107      REAL(wp) ::   zsto, zjulian, zout,   &  ! temporary scalars
108         &          zindh, zinda, zindb, ztmu
109      REAL(wp), DIMENSION(1)                ::   zdept
110      REAL(wp), POINTER, DIMENSION(:,:)     ::   zfield
111      !!-------------------------------------------------------------------
112
113      CALL wrk_alloc( jpi, jpj, zfield )
114                                                 !--------------------!
115      IF( kt == nit000 ) THEN                    !   Initialisation   !
116         !                                       !--------------------!
117
118         CALL lim_wri_init_2 
119                           
120         zsto     = rdt_ice
121         IF( ln_mskland )   THEN   ;   clop = "ave(only(x))"   ! put 1.e+20 on land (very expensive!!)
122         ELSE                      ;   clop = "ave(x)"         ! no use of the mask value (require less cpu time)
123         ENDIF
124         zout     = nwrite * rdt_ice / nn_fsbc
125         niter    = ( nit000 - 1 ) / nn_fsbc
126         zdept(1) = 0.
127         
128         CALL ymds2ju ( nyear, nmonth, nday, rdt, zjulian )
129         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
130         CALL dia_nam ( clhstnam, nwrite, 'icemod' )
131         CALL histbeg ( clhstnam, jpi, glamt, jpj, gphit,    &
132            &           1, jpi, 1, jpj, niter, zjulian, rdt_ice, nhorid, nice , domain_id=nidom, snc4chunks=snc4set)
133         CALL histvert( nice, "deptht", "Vertical T levels", "m", 1, zdept, ndepid, "down")
134         CALL wheneq  ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim)
135         
136         DO jf = 1, noumef
137            IF( nc(jf) == 1 )   CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj   &
138               &                                  , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout )
139         END DO
140         CALL histend( nice, snc4set )
141         !
142      ENDIF
143      !                                          !--------------------!
144      !                                          !   Cumulate at kt   !
145      !                                          !--------------------!
146
147      !-- Store instantaneous values in zcmo
148     
149      zcmo(:,:, 1:jpnoumax ) = 0.e0 
150      DO jj = 2 , jpjm1
151         DO ji = 1 , jpim1   ! NO vector opt.
152            zindh  = MAX( zzero , SIGN( zone , hicif(ji,jj) * (1.0 - frld(ji,jj) ) - 0.10 ) )
153            zinda  = MAX( zzero , SIGN( zone , ( 1.0 - frld(ji,jj) ) - 0.10 ) )
154            zindb  = zindh * zinda
155            ztmu   = MAX( 0.5 * zone , ( tmu(ji,jj) + tmu(ji+1,jj) + tmu(ji,jj+1) + tmu(ji+1,jj+1) ) ) 
156            zcmo(ji,jj,1)  = hsnif (ji,jj)
157            zcmo(ji,jj,2)  = hicif (ji,jj)
158            zcmo(ji,jj,3)  = hicifp(ji,jj)
159            zcmo(ji,jj,4)  = frld  (ji,jj)
160            zcmo(ji,jj,5)  = sist  (ji,jj)
161            zcmo(ji,jj,6)  = fbif  (ji,jj)
162           IF (lk_lim2_vp) THEN
163            zcmo(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * tmu(ji,jj  ) + u_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
164                                      + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
165                                  / ztmu 
166
167            zcmo(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * tmu(ji,jj  ) + v_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
168                                      + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
169                                  / ztmu
170           ELSE
171
172            zcmo(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * tmu(ji,jj)                       &
173             &                        + u_ice(ji-1,jj) * tmu(ji-1,jj) )                   &
174             &                    / 2.0
175            zcmo(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * tmv(ji,jj)                       &
176             &                        + v_ice(ji,jj-1) * tmv(ji,jj-1) )                   &
177             &                    / 2.0
178
179           ENDIF
180            zcmo(ji,jj,9)  = sst_m(ji,jj)
181            zcmo(ji,jj,10) = sss_m(ji,jj)
182            zcmo(ji,jj,11) = qns(ji,jj) + qsr(ji,jj)
183            zcmo(ji,jj,12) = qsr(ji,jj)
184            zcmo(ji,jj,13) = qns(ji,jj)
185            ! See thersf for the coefficient
186            zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce    !!gm ???
187            zcmo(ji,jj,15) = utau_ice(ji,jj)
188            zcmo(ji,jj,16) = vtau_ice(ji,jj)
189            zcmo(ji,jj,17) = qsr_ice(ji,jj,1)
190            zcmo(ji,jj,18) = qns_ice(ji,jj,1)
191            zcmo(ji,jj,19) = sprecip(ji,jj)
192         END DO
193      END DO
194      !
195      ! Write the netcdf file
196      !
197      niter = niter + 1
198      DO jf = 1 , noumef
199         zfield(:,:) = zcmo(:,:,jf) * cmulti(jf) + cadd(jf) * tmask(:,:,1)
200         SELECT CASE ( jf )
201         CASE ( 7, 8, 15, 16, 20, 21 )  ! velocity or stress fields (vectors)
202            CALL lbc_lnk( zfield, 'T', -1. )
203         CASE DEFAULT                   ! scalar fields
204            CALL lbc_lnk( zfield, 'T',  1. )
205         END SELECT
206
207         IF( nc(jf) == 1 )   CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 )
208
209      END DO
210
211      IF( ( nn_fsbc * niter ) >= nitend )   CALL histclo( nice ) 
212
213      CALL wrk_dealloc( jpi, jpj, zfield )
214      !
215   END SUBROUTINE lim_wri_2
216     
217#endif     
218
219   SUBROUTINE lim_wri_init_2
220      !!-------------------------------------------------------------------
221      !!                    ***   ROUTINE lim_wri_init_2  ***
222      !!               
223      !! ** Purpose :   intialisation of LIM sea-ice output
224      !!
225      !! ** Method  : Read the namicewri namelist and check the parameter
226      !!       values called at the first timestep (nit000)
227      !!
228      !! ** input   :   Namelist namicewri
229      !!-------------------------------------------------------------------
230      INTEGER ::   nf      ! ???
231      TYPE FIELD 
232         CHARACTER(len = 35) :: ztitle 
233         CHARACTER(len = 8 ) :: zname         
234         CHARACTER(len = 8 ) :: zunit
235         INTEGER             :: znc   
236         REAL                :: zcmulti 
237         REAL                :: zcadd       
238      END TYPE FIELD
239      TYPE(FIELD) ::  &
240         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,   &
241         field_7 , field_8 , field_9 , field_10, field_11, field_12,   &
242         field_13, field_14, field_15, field_16, field_17, field_18,   &
243         field_19
244      TYPE(FIELD) , DIMENSION(jpnoumax) :: zfield
245
246      NAMELIST/namiceout/ noumef, &
247         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,   &
248         field_7 , field_8 , field_9 , field_10, field_11, field_12,   &
249         field_13, field_14, field_15, field_16, field_17, field_18,   &
250         field_19
251      !!-------------------------------------------------------------------
252      !
253      IF( lim_wri_alloc_2() /= 0 ) THEN      ! allocate lim_wri arrrays
254         CALL ctl_stop( 'STOP', 'lim_wri_init_2 : unable to allocate standard arrays' )   ;   RETURN
255      ENDIF
256
257      REWIND ( numnam_ice )                ! Read Namelist namicewri
258      READ   ( numnam_ice  , namiceout )
259     
260      zfield( 1) = field_1
261      zfield( 2) = field_2
262      zfield( 3) = field_3
263      zfield( 4) = field_4
264      zfield( 5) = field_5
265      zfield( 6) = field_6
266      zfield( 7) = field_7
267      zfield( 8) = field_8
268      zfield( 9) = field_9
269      zfield(10) = field_10
270      zfield(11) = field_11
271      zfield(12) = field_12
272      zfield(13) = field_13
273      zfield(14) = field_14
274      zfield(15) = field_15
275      zfield(16) = field_16
276      zfield(17) = field_17
277      zfield(18) = field_18
278      zfield(19) = field_19
279     
280      DO nf = 1, noumef
281         titn  (nf) = zfield(nf)%ztitle
282         nam   (nf) = zfield(nf)%zname
283         uni   (nf) = zfield(nf)%zunit
284         nc    (nf) = zfield(nf)%znc
285         cmulti(nf) = zfield(nf)%zcmulti
286         cadd  (nf) = zfield(nf)%zcadd
287      END DO
288
289      IF(lwp) THEN
290         WRITE(numout,*)
291         WRITE(numout,*) 'lim_wri_init_2 : Ice parameters for outputs'
292         WRITE(numout,*) '~~~~~~~~~~~~~~'
293         WRITE(numout,*) '    number of fields to be stored         noumef = ', noumef
294         WRITE(numout,*) '           title                            name     unit      Saving (1/0) ',   &
295            &            '    multiplicative constant       additive constant '
296         DO nf = 1 , noumef         
297            WRITE(numout,*) '   ', titn(nf), '   ', nam(nf),'      ', uni(nf),'  ', nc(nf),'        ', cmulti(nf),   &
298               &       '        ', cadd(nf)
299         END DO
300      ENDIF
301      !   
302   END SUBROUTINE lim_wri_init_2
303
304#endif
305
306   SUBROUTINE lim_wri_state_2( kt, kid, kh_i )
307      !!---------------------------------------------------------------------
308      !!                 ***  ROUTINE lim_wri_state_2  ***
309      !!       
310      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
311      !!      the instantaneous ice state and forcing fields for ice model
312      !!        Used to find errors in the initial state or save the last
313      !!      ocean state in case of abnormal end of a simulation
314      !!
315      !! History :
316      !!   2.0  !  2009-06  (B. Lemaire)
317      !!----------------------------------------------------------------------
318      INTEGER, INTENT( in ) ::   kt               ! ocean time-step index)
319      INTEGER, INTENT( in ) ::   kid , kh_i       
320      !!----------------------------------------------------------------------
321
322      CALL histdef( kid, "isnowthi", "Snow thickness"          , "m"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
323      CALL histdef( kid, "iicethic", "Ice thickness"           , "m"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
324      CALL histdef( kid, "iiceprod", "Ice produced"            , "m/kt"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
325      CALL histdef( kid, "ileadfra", "Ice concentration"       , "-"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
326      CALL histdef( kid, "iicetemp", "Ice temperature"         , "K"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
327      CALL histdef( kid, "ioceflxb", "flux at ice base"        , "w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
328      CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)"   , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
329      CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)"   , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
330      CALL histdef( kid, "isstempe", "Sea surface temperature" , "C"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
331      CALL histdef( kid, "isssalin", "Sea surface salinity"    , "PSU"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
332      CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
333      CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
334      CALL histdef( kid, "iicesflx", "Solar flux over ice"     , "w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
335      CALL histdef( kid, "iicenflx", "Non-solar flux over ice" , "w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
336      CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
337
338      CALL histend( kid, snc4set )   ! end of the file definition
339
340      CALL histwrite( kid, "isnowthi", kt, hsnif          , jpi*jpj, (/1/) )   
341      CALL histwrite( kid, "iicethic", kt, hicif          , jpi*jpj, (/1/) )   
342      CALL histwrite( kid, "iiceprod", kt, hicifp         , jpi*jpj, (/1/) )   
343      CALL histwrite( kid, "ileadfra", kt, 1. - frld(:,:) , jpi*jpj, (/1/) )
344      CALL histwrite( kid, "iicetemp", kt, sist(:,:) - rt0, jpi*jpj, (/1/) )
345      CALL histwrite( kid, "ioceflxb", kt, fbif           , jpi*jpj, (/1/) )
346      CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) )
347      CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) )
348      CALL histwrite( kid, "isstempe", kt, sst_m          , jpi*jpj, (/1/) )
349      CALL histwrite( kid, "isssalin", kt, sss_m          , jpi*jpj, (/1/) )
350      CALL histwrite( kid, "iicestru", kt, utau_ice       , jpi*jpj, (/1/) )
351      CALL histwrite( kid, "iicestrv", kt, vtau_ice       , jpi*jpj, (/1/) )
352      CALL histwrite( kid, "iicesflx", kt, qsr_ice(:,:,1) , jpi*jpj, (/1/) )
353      CALL histwrite( kid, "iicenflx", kt, qns_ice(:,:,1) , jpi*jpj, (/1/) )
354      CALL histwrite( kid, "isnowpre", kt, sprecip        , jpi*jpj, (/1/) )
355
356    END SUBROUTINE lim_wri_state_2
357
358#else
359   !!----------------------------------------------------------------------
360   !!   Default option :         Empty module      NO LIM 2.0 sea-ice model
361   !!----------------------------------------------------------------------
362CONTAINS
363   SUBROUTINE lim_wri_2          ! Empty routine
364   END SUBROUTINE lim_wri_2
365#endif
366
367   !!======================================================================
368END MODULE limwri_2
Note: See TracBrowser for help on using the repository browser.