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.
limdia.F90 on Ticket #741 – Attachment – NEMO

Ticket #741: limdia.F90

File limdia.F90, 29.7 KB (added by vancop, 14 years ago)

modified routine

Line 
1MODULE limdia
2   !!======================================================================
3   !!                       ***  MODULE limdia   ***
4   !!                      diagnostics of ice model
5   !!======================================================================
6   !! To add a new field
7   !! 1) in lim_dia : add its definition for both hemispheres if wished
8   !! 2) add the new titles in lim_dia_init
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3' :                                   LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_dia      : computation of the time evolution of keys var.
15   !!   lim_dia_init : initialization and namelist read
16   !!----------------------------------------------------------------------
17   !! * Modules used
18   USE phycst
19   USE in_out_manager
20   USE par_ice         ! ice parameters
21   USE ice_oce         ! ice variables
22   USE daymod
23   USE dom_ice
24   USE ice
25   USE iceini
26   USE limistate
27   USE dom_oce
28   USE sbc_oce         ! Surface boundary condition: ocean fields
29   USE dom_oce
30   USE lib_mpp
31   IMPLICIT NONE
32   PRIVATE
33
34   !! * Routine accessibility
35   PUBLIC lim_dia       ! called by ice_step
36
37   !! * Shared module variables
38   INTEGER, PUBLIC  ::  &
39      ntmoy   = 1 ,     &  !: instantaneous values of ice evolution or averaging ntmoy
40      ninfo   = 1          !: frequency of ouputs on file ice_evolu in case of averaging
41
42   !! * Module variables
43   INTEGER, PARAMETER ::   &  ! Parameters for outputs to files "evolu"
44      jpinfmx = 100         ,    &  ! maximum number of key variables
45      jpchinf = 5           ,    &  ! ???
46      jpchsep = jpchinf + 2         ! ???
47
48   INTEGER          ::  &
49      nfrinf  = 4    ,  &  ! number of variables written in one line
50      nferme ,          &  ! last time step at which the var. are written on file
51      nvinfo ,          &  ! number of total variables
52      nbvt   ,          &  ! number of time variables
53      naveg                ! number of step for accumulation before averaging
54
55   CHARACTER(len=8) ::  &
56      fmtinf  = '1PE13.5 ' ! format of the output values 
57
58   CHARACTER(len=30) :: &
59      fmtw  ,           &  ! formats
60      fmtr  ,           &  ! ???
61      fmtitr               ! ???
62
63   CHARACTER(len=jpchsep), DIMENSION(jpinfmx) ::   &
64      titvar               ! title of key variables
65
66   REAL(wp), DIMENSION(jpinfmx) ::  &
67      vinfom               ! temporary working space
68   REAL(wp) :: &
69      epsi06 = 1.e-06
70   REAL(wp), DIMENSION(jpi,jpj) ::   &
71      aire                 ! masked grid cell area
72
73   !! * Substitutions
74#  include "vectopt_loop_substitute.h90"
75   !!----------------------------------------------------------------------
76   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)
77   !! $Id: limdia.F90 1156 2008-06-26 16:06:45Z rblod $
78   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
79   !!----------------------------------------------------------------------
80
81CONTAINS
82
83   SUBROUTINE lim_dia
84      !!--------------------------------------------------------------------
85      !!                  ***  ROUTINE lim_dia  ***
86      !!   
87      !! ** Purpose : Computation and outputs on file ice.evolu
88      !!      the temporal evolution of some key variables
89      !!
90      !! History :
91      !!   8.0  !  97-06  (Louvain-La-Neuve)  Original code
92      !!   8.5  !  02-09  (C. Ethe , G. Madec )  F90: Free form and module
93      !!   9.0  !  05-06  (M. Vancoppenolle LLN) Rehabilitation and fixing a few bugs
94      !!           01-07  (M. Vancoppenolle LLN) added new fields + fram strait
95      !!                                         export
96      !!-------------------------------------------------------------------
97      !! * Local variables
98      INTEGER  ::   jv,ji,jj,jl ! dummy loop indices
99      REAL(wp), DIMENSION(jpinfmx) ::  & 
100         vinfor           ! temporary working space
101      REAL(wp) ::    &
102         zshift_date   , & ! date from the minimum ice extent
103         zday, zday_min, & ! current day, day of minimum extent
104         zafy, zamy,     & ! temporary area of fy and my ice
105         zindb
106      !!-------------------------------------------------------------------
107
108      ! 0) date from the minimum of ice extent
109      !---------------------------------------
110      zday_min = 273.0        ! zday_min = date of minimum extent, here September 30th
111      zday = FLOAT(numit-nit000) * rdt_ice / ( 86400.0 * FLOAT(nn_fsbc) )
112      IF (zday.GT.zday_min) THEN
113         zshift_date  =  zday - zday_min
114      ELSE
115         zshift_date  =  zday - (365.0 - zday_min)
116      ENDIF
117
118      IF( numit == nstart )   CALL lim_dia_init   ! initialisation of ice_evolu file     
119
120      ! temporal diagnostics
121      vinfor(1) = REAL(numit)
122      vinfor(2) = nyear
123
124      ! put everything to zero
125      DO jv = nbvt + 1, nvinfo
126         vinfor(jv) = 0.0
127      END DO
128
129      !!-------------------------------------------------------------------
130      !! 1) Northern hemisphere
131      !!-------------------------------------------------------------------
132      !! 1.1) Diagnostics independent on age
133      !!------------------------------------
134      DO jj = njeq, jpjm1
135         DO ji = fs_2, fs_jpim1   ! vector opt.
136            IF( tms(ji,jj) == 1 ) THEN
137               vinfor(3)  = vinfor(3)  + at_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice area
138               IF (at_i(ji,jj).GT.0.15) vinfor(5) = vinfor(5) + aire(ji,jj) / 1.0e12 !ice extent
139               vinfor(7)  = vinfor(7)  + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume
140               vinfor(9)  = vinfor(9)  + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume
141               vinfor(15) = vinfor(15) + ot_i(ji,jj) *vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age
142               vinfor(29) = vinfor(29) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity
143               ! the computation of this diagnostic is not reliable
144               vinfor(31) = vinfor(31) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 
145                  v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 
146               vinfor(53) = vinfor(53) + emps(ji,jj)*aire(ji,jj) / 1.0e12 !salt flux
147               vinfor(55) = vinfor(55) + fsbri(ji,jj)*aire(ji,jj) / 1.0e12 !brine drainage flux
148               vinfor(57) = vinfor(57) + fseqv(ji,jj)*aire(ji,jj) / 1.0e12 !equivalent salt flux
149               vinfor(59) = vinfor(59) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SST
150               vinfor(61) = vinfor(61) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SSS
151               vinfor(65) = vinfor(65) + et_s(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12  ! snow temperature
152               vinfor(67) = vinfor(67) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12       ! ice heat content
153               vinfor(69) = vinfor(69) + v_i(ji,jj,1)*aire(ji,jj) / 1.0e12 !ice volume
154               vinfor(71) = vinfor(71) + v_i(ji,jj,2)*aire(ji,jj) / 1.0e12 !ice volume
155               vinfor(73) = vinfor(73) + v_i(ji,jj,3)*aire(ji,jj) / 1.0e12 !ice volume
156               vinfor(75) = vinfor(75) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume
157               vinfor(77) = vinfor(77) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume
158               vinfor(79) = 0.0
159               vinfor(81) = vinfor(81) + emp(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux
160            ENDIF
161         END DO
162      END DO
163
164      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
165         DO jj = njeq, jpjm1
166            DO ji = fs_2, fs_jpim1   ! vector opt.
167               IF( tms(ji,jj) == 1 ) THEN
168                  vinfor(11) = vinfor(11) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !undef def ice volume
169               ENDIF
170            END DO
171         END DO
172      END DO
173
174      vinfor(13) = 0.0
175
176      vinfor(15) = vinfor(15) / MAX(vinfor(7),epsi06) ! these have to be divided by total ice volume to have the
177      vinfor(29) = vinfor(29) / MAX(vinfor(7),epsi06) ! right value
178      vinfor(31) = SQRT( vinfor(31) / MAX( vinfor(7) , epsi06 ) )
179      vinfor(67) = vinfor(67) / MAX(vinfor(7),epsi06)
180
181      vinfor(53) = vinfor(53) / MAX(vinfor(5),epsi06) ! these have to be divided by total ice extent to have the
182      vinfor(55) = vinfor(55) / MAX(vinfor(5),epsi06) ! right value
183      vinfor(57) = vinfor(57) / MAX(vinfor(5),epsi06) !
184      vinfor(79) = vinfor(79) / MAX(vinfor(5),epsi06) !
185
186      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(3))) !
187      vinfor(59) = zindb*vinfor(59) / MAX(vinfor(3),epsi06) ! divide by ice area
188      vinfor(61) = zindb*vinfor(61) / MAX(vinfor(3),epsi06) !
189
190      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(9))) !
191      vinfor(65) = zindb*vinfor(65) / MAX(vinfor(9),epsi06) ! divide it by snow volume
192
193
194      DO jl = 1, jpl
195         DO jj = njeq, jpjm1
196            DO ji = fs_2, fs_jpim1   ! vector opt.
197               IF( tms(ji,jj) == 1 ) THEN
198                  vinfor(33) = vinfor(33) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume
199                  vinfor(35) = vinfor(35) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume
200               ENDIF
201            END DO
202         END DO
203      END DO
204
205      DO jj = njeq, jpjm1
206         DO ji = fs_2, fs_jpim1   ! vector opt.
207            IF( tms(ji,jj) == 1 ) THEN
208               vinfor(37) = vinfor(37) + diag_sni_gr(ji,jj)*aire(ji,jj) / 1.0e12 !th growth rates
209               vinfor(39) = vinfor(39) + diag_lat_gr(ji,jj)*aire(ji,jj) / 1.0e12 
210               vinfor(41) = vinfor(41) + diag_bot_gr(ji,jj)*aire(ji,jj) / 1.0e12
211               vinfor(43) = vinfor(43) + diag_dyn_gr(ji,jj)*aire(ji,jj) / 1.0e12 
212               vinfor(45) = vinfor(45) + dv_dt_thd(ji,jj,5)*aire(ji,jj) / 1.0e12
213               vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) / 1.0e12 / rdt_ice ! volume acc in OW
214            ENDIF
215         END DO
216      END DO
217
218      DO jl = 1, jpl
219         DO jj = njeq, jpjm1
220            DO ji = fs_2, fs_jpim1   ! vector opt.
221               IF( tms(ji,jj) == 1 ) THEN
222                  vinfor(63) = vinfor(63) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12
223               ENDIF
224            END DO
225         END DO
226      END DO
227      vinfor(63) = vinfor(63) / MAX(vinfor(3),epsi06) ! these have to be divided by total ice area
228
229      !! 1.2) Diagnostics dependent on age
230      !!------------------------------------
231      DO jj = njeq, jpjm1
232         DO ji = fs_2, fs_jpim1   ! vector opt.
233            IF( tms(ji,jj) == 1 ) THEN
234               zafy = 0.0
235               zamy = 0.0
236               DO jl = 1, jpl
237                  IF ((o_i(ji,jj,jl) - zshift_date).LT.0.0) THEN
238                     vinfor(17) = vinfor(17) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice area
239                     vinfor(25) = vinfor(25) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice volume
240                     vinfor(49) = vinfor(49) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity
241                     zafy = zafy + a_i(ji,jj,jl)
242                  ENDIF
243                  IF ((o_i(ji,jj,jl) - zshift_date).GT.0.0) THEN
244                     vinfor(19) = vinfor(19) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12    ! MY ice area
245                     vinfor(27) = vinfor(27) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! MY ice volume
246                     vinfor(51) = vinfor(51) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !MY ice salinity
247                     zamy = zamy + a_i(ji,jj,jl)
248                  ENDIF
249               END DO
250               IF ((at_i(ji,jj).GT.0.15).AND.(zafy.GT.zamy)) THEN
251                  vinfor(21) = vinfor(21) + aire(ji,jj) / 1.0e12 ! Seasonal ice extent
252               ENDIF
253               IF ((at_i(ji,jj).GT.0.15).AND.(zafy.LE.zamy)) THEN
254                  vinfor(23) = vinfor(23) + aire(ji,jj) / 1.0e12 ! Perennial ice extent
255               ENDIF
256            ENDIF
257         END DO
258      END DO
259      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(25))) !=0 if no multiyear ice 1 if yes
260      vinfor(49) = zindb*vinfor(49) / MAX(vinfor(25),epsi06)
261      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(27))) !=0 if no multiyear ice 1 if yes
262      vinfor(51) = zindb*vinfor(51) / MAX(vinfor(27),epsi06)
263
264      !! Fram Strait Export
265      !! 83 = area export
266      !! 84 = volume export
267      !! Fram strait in ORCA2 = 5 points
268      !! export = -v_ice*e1t*ddtb*at_i or -v_ice*e1t*ddtb*at_i*h_i
269      jj = 136 ! C grid
270      vinfor(83) = 0.0
271      vinfor(84) = 0.0
272      DO ji = 134, 138
273         vinfor(83) = vinfor(83) - v_ice(ji,jj) * & 
274            e1t(ji,jj)*at_i(ji,jj)*rdt_ice / 1.0e12
275         vinfor(84) = vinfor(84) - v_ice(ji,jj) * & 
276            e1t(ji,jj)*vt_i(ji,jj)*rdt_ice / 1.0e12
277      END DO
278
279      !!-------------------------------------------------------------------
280      !! 2) Southern hemisphere
281      !!-------------------------------------------------------------------
282      !! 2.1) Diagnostics independent on age
283      !!------------------------------------
284      DO jj = 2, njeqm1
285         DO ji = fs_2, fs_jpim1   ! vector opt.
286            IF( tms(ji,jj) == 1 ) THEN
287               vinfor(4)  = vinfor(4)  + at_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice area
288               IF (at_i(ji,jj).GT.0.15) vinfor(6) = vinfor(6) + aire(ji,jj) / 1.0e12 !ice extent
289               vinfor(8)  = vinfor(8)  + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume
290               vinfor(10) = vinfor(10) + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume
291               vinfor(16) = vinfor(16) + ot_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age
292               vinfor(30) = vinfor(30) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity
293               ! this diagnostic is not well computed (weighted by vol instead
294               ! of area)
295               vinfor(32) = vinfor(32) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 
296                  v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel
297               vinfor(54) = vinfor(54) + at_i(ji,jj)*emps(ji,jj)*aire(ji,jj) / 1.0e12 ! Total salt flux
298               vinfor(56) = vinfor(56) + at_i(ji,jj)*fsbri(ji,jj)*aire(ji,jj) / 1.0e12 ! Brine drainage salt flux
299               vinfor(58) = vinfor(58) + at_i(ji,jj)*fseqv(ji,jj)*aire(ji,jj) / 1.0e12 ! Equivalent salt flux
300               vinfor(60) = vinfor(60) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SST
301               vinfor(62) = vinfor(62) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SSS
302               vinfor(66) = vinfor(66) + et_s(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! snow temperature
303               vinfor(68) = vinfor(68) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! ice enthalpy
304               vinfor(70) = vinfor(70) + v_i(ji,jj,1)*aire(ji,jj) / 1.0e12 !ice volume
305               vinfor(72) = vinfor(72) + v_i(ji,jj,2)*aire(ji,jj) / 1.0e12 !ice volume
306               vinfor(74) = vinfor(74) + v_i(ji,jj,3)*aire(ji,jj) / 1.0e12 !ice volume
307               vinfor(76) = vinfor(76) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume
308               vinfor(78) = vinfor(78) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume
309               vinfor(80) = 0.0
310               vinfor(82) = vinfor(82) + emp(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux
311            ENDIF
312         END DO
313      END DO
314
315      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
316         DO jj = 2, njeqm1
317            DO ji = fs_2, fs_jpim1   ! vector opt.
318               vinfor(12) = vinfor(12) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !undef def ice volume
319            END DO
320         END DO
321      END DO
322
323      vinfor(14) = 0.0
324
325      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(8))) 
326      vinfor(16) = zindb * vinfor(16) / MAX(vinfor(8),epsi06) ! these have to be divided by ice vol
327      vinfor(30) = zindb * vinfor(30) / MAX(vinfor(8),epsi06) !
328      vinfor(32) = zindb * SQRT( vinfor(32) / MAX( vinfor(8) , epsi06 ) )
329      vinfor(68) = zindb * vinfor(68) / MAX(vinfor(8),epsi06) !
330
331      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(6))) 
332      vinfor(54) = zindb * vinfor(54) / MAX(vinfor(6),epsi06) ! these have to be divided by ice extt
333      vinfor(56) = zindb * vinfor(56) / MAX(vinfor(6),epsi06) !
334      vinfor(58) = zindb * vinfor(58) / MAX(vinfor(6),epsi06) !
335      vinfor(80) = zindb * vinfor(80) / MAX(vinfor(6),epsi06) !
336      !      vinfor(84) = vinfor(84) / vinfor(6) !
337
338      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4))) !
339      vinfor(60) = zindb*vinfor(60) / ( MAX(vinfor(4), epsi06) ) ! divide by ice area
340      vinfor(62) = zindb*vinfor(62) / ( MAX(vinfor(4), epsi06) ) !
341
342      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(10))) !
343      vinfor(66) = zindb*vinfor(66) / MAX(vinfor(10),epsi06) ! divide it by snow volume
344
345      DO jl = 1, jpl
346         DO jj = 2, njeqm1
347            DO ji = fs_2, fs_jpim1   ! vector opt.
348               IF( tms(ji,jj) == 1 ) THEN
349                  vinfor(34) = vinfor(34) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume
350                  vinfor(36) = vinfor(36) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume
351               ENDIF
352            END DO
353         END DO
354      END DO
355
356      DO jj = 2, njeqm1
357         DO ji = fs_2, fs_jpim1   ! vector opt.
358            IF( tms(ji,jj) == 1 ) THEN
359               vinfor(38) = vinfor(38) + diag_sni_gr(ji,jj)*aire(ji,jj) / 1.0e12 !th growth rates
360               vinfor(40) = vinfor(40) + diag_lat_gr(ji,jj)*aire(ji,jj) / 1.0e12 
361               vinfor(42) = vinfor(42) + diag_bot_gr(ji,jj)*aire(ji,jj) / 1.0e12
362               vinfor(44) = vinfor(44) + diag_dyn_gr(ji,jj)*aire(ji,jj) / 1.0e12 
363               vinfor(46) = vinfor(46) + dv_dt_thd(ji,jj,5)*aire(ji,jj) / 1.0e12
364               vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) / 1.0e12 / rdt_ice ! volume acc in OW
365            ENDIF
366         END DO
367      END DO
368
369
370      DO jl = 1, jpl
371         DO jj = 2, njeqm1
372            DO ji = fs_2, fs_jpim1   ! vector opt.
373               IF( tms(ji,jj) == 1 ) THEN
374                  vinfor(64) = vinfor(64) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12
375               ENDIF
376            END DO
377         END DO
378      END DO
379      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4))) !
380      vinfor(64) = zindb * vinfor(64) / MAX(vinfor(4),epsi06) ! divide by ice extt
381      !! 2.2) Diagnostics dependent on age
382      !!------------------------------------
383      DO jj = 2, njeqm1
384         DO ji = fs_2, fs_jpim1   ! vector opt.
385            IF( tms(ji,jj) == 1 ) THEN
386               zafy = 0.0
387               zamy = 0.0
388               DO jl = 1, jpl
389                  IF ((o_i(ji,jj,jl) - zshift_date).LT.0.0) THEN
390                     vinfor(18) = vinfor(18) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice area
391                     vinfor(26) = vinfor(26) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice volume
392                     zafy = zafy + a_i(ji,jj,jl)
393                     vinfor(50) = vinfor(50) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity
394                  ENDIF
395                  IF ((o_i(ji,jj,jl) - zshift_date).GT.0.0) THEN
396                     vinfor(20) = vinfor(20) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12    ! MY ice area
397                     vinfor(28) = vinfor(28) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12
398                     vinfor(52) = vinfor(52) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity
399                     zamy = zamy + a_i(ji,jj,jl)
400                  ENDIF
401               END DO ! jl
402               IF ((at_i(ji,jj).GT.0.15).AND.(zafy.GT.zamy)) THEN
403                  vinfor(22) = vinfor(22) + aire(ji,jj) / 1.0e12 ! Seasonal ice extent
404               ENDIF
405               IF ((at_i(ji,jj).GT.0.15).AND.(zafy.LE.zamy)) THEN
406                  vinfor(24) = vinfor(24) + aire(ji,jj) / 1.0e12 ! Perennial ice extent
407               ENDIF
408            ENDIF ! tms
409         END DO ! jj
410      END DO ! ji
411      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(26))) !=0 if no multiyear ice 1 if yes
412      vinfor(50) = zindb*vinfor(50) / MAX(vinfor(26),epsi06)
413      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(28))) !=0 if no multiyear ice 1 if yes
414      vinfor(52) = zindb*vinfor(52) / MAX(vinfor(28),epsi06)
415
416      !  Accumulation before averaging
417      DO jv = 1, nvinfo
418         vinfom(jv) = vinfom(jv) + vinfor(jv)
419      END DO
420      naveg = naveg + 1 
421
422      ! oututs on file ice_evolu   
423      !MV      IF( MOD( numit , ninfo ) == 0 ) THEN
424      WRITE(numevo_ice,fmtw) ( titvar(jv), vinfom(jv)/naveg, jv = 1, nvinfo )
425      naveg = 0
426      DO jv = 1, nvinfo
427         vinfom(jv)=0.0
428      END DO
429      !MV      ENDIF
430
431   END SUBROUTINE lim_dia
432
433   SUBROUTINE lim_dia_init
434      !!-------------------------------------------------------------------
435      !!                  ***  ROUTINE lim_dia_init  ***
436      !!             
437      !! ** Purpose : Preparation of the file ice_evolu for the output of
438      !!      the temporal evolution of key variables
439      !!
440      !! ** input   : Namelist namicedia
441      !!
442      !! history :
443      !!  8.5  ! 03-08 (C. Ethe) original code
444      !!  9.0  ! 08-03 (M. Vancoppenolle) LIM3
445      !!-------------------------------------------------------------------
446      NAMELIST/namicedia/fmtinf, nfrinf, ninfo, ntmoy
447
448      INTEGER  ::   jv   ,     &  ! dummy loop indice
449         &          ntot ,     &
450         &          ndeb ,     &
451         &          irecl
452
453      REAL(wp) ::   zxx0, zxx1    ! temporary scalars
454
455      CHARACTER(len=jpchinf) ::   titinf
456      CHARACTER(len=50)      ::   clname
457      !!-------------------------------------------------------------------
458
459
460      ! Read Namelist namicedia
461      REWIND ( numnam_ice )
462      READ   ( numnam_ice  , namicedia )
463      IF(lwp) THEN
464         WRITE(numout,*)
465         WRITE(numout,*) 'lim_dia_init : ice parameters for ice diagnostics '
466         WRITE(numout,*) '~~~~~~~~~~~~'
467         WRITE(numout,*) '   format of the output values                                 fmtinf = ', fmtinf
468         WRITE(numout,*) '   number of variables written in one line                     nfrinf = ', nfrinf 
469         WRITE(numout,*) '   Instantaneous values of ice evolution or averaging          ntmoy  = ', ntmoy
470         WRITE(numout,*) '   frequency of ouputs on file ice_evolu in case of averaging  ninfo  = ', ninfo
471      ENDIF
472
473      ! masked grid cell area
474      aire(:,:) = area(:,:) * tms(:,:)
475
476      ! Titles of ice key variables :
477      titvar(1) = 'NoIt'  ! iteration number
478      titvar(2) = 'T yr'  ! time step in years
479      nbvt = 2            ! number of time variables
480
481      titvar(3) = 'AI_N'  ! sea ice area in the northern Hemisp.(10^12 km2)
482      titvar(4) = 'AI_S'  ! sea ice area in the southern Hemisp.(10^12 km2)
483      titvar(5) = 'EI_N'  ! sea ice extent (15%) in the northern Hemisp.(10^12 km2)
484      titvar(6) = 'EI_S'  ! sea ice extent (15%) in the southern Hemisp.(10^12 km2)
485      titvar(7) = 'VI_N'  ! sea ice volume in the northern Hemisp.(10^3 km3)
486      titvar(8) = 'VI_S'  ! sea ice volume in the southern Hemisp.(10^3 km3)
487      titvar(9) = 'VS_N'  ! snow volume over sea ice in the northern Hemisp.(10^3 km3)
488      titvar(10)= 'VS_S'  ! snow volume over sea ice in the northern Hemisp.(10^3 km3)
489      titvar(11)= 'VuIN'  ! undeformed sea ice volume in the northern Hemisp.(10^3 km3)
490      titvar(12)= 'VuIS'  ! undeformed sea ice volume in the southern Hemisp.(10^3 km3)
491      titvar(13)= 'VdIN'  ! deformed sea ice volume in the northern Hemisp.(10^3 km3)
492      titvar(14)= 'VdIS'  ! deformed sea ice volume in the southern Hemisp.(10^3 km3)
493      titvar(15)= 'OI_N'  ! sea ice mean age in the northern Hemisp.(years)
494      titvar(16)= 'OI_S'  ! sea ice mean age in the southern Hemisp.(years)
495      titvar(17)= 'AFYN'  ! total FY ice area northern Hemisp.(10^12 km2)
496      titvar(18)= 'AFYS'  ! total FY ice area southern Hemisp.(10^12 km2)
497      titvar(19)= 'AMYN'  ! total MY ice area northern Hemisp.(10^12 km2)
498      titvar(20)= 'AMYS'  ! total MY ice area southern Hemisp.(10^12 km2)
499      titvar(21)= 'EFYN'  ! total FY ice extent northern Hemisp.(10^12 km2) (with more 50% FY ice)
500      titvar(22)= 'EFYS'  ! total FY ice extent southern Hemisp.(10^12 km2) (with more 50% FY ice)
501      titvar(23)= 'EMYN'  ! total MY ice extent northern Hemisp.(10^12 km2) (with more 50% MY ice)
502      titvar(24)= 'EMYS'  ! total MY ice extent southern Hemisp.(10^12 km2) (with more 50% MY ice)
503      titvar(25)= 'VFYN'  ! total undeformed FY ice volume northern Hemisp.(10^3 km3)
504      titvar(26)= 'VFYS'  ! total undeformed FY ice volume southern Hemisp.(10^3 km3)
505      titvar(27)= 'VMYN'  ! total undeformed MY ice volume northern Hemisp.(10^3 km3)
506      titvar(28)= 'VMYS'  ! total undeformed MY ice volume southern Hemisp.(10^3 km3)
507      titvar(29)= 'IS_N'  ! sea ice mean salinity in the northern hemisphere (ppt) 
508      titvar(30)= 'IS_S'  ! sea ice mean salinity in the southern hemisphere (ppt) 
509      titvar(31)= 'IVeN'  ! sea ice mean velocity in the northern hemisphere (m/s)
510      titvar(32)= 'IVeS'  ! sea ice mean velocity in the southern hemisphere (m/s)
511      titvar(33)= 'DVDN'  ! variation of sea ice volume due to dynamics in the northern hemisphere
512      titvar(34)= 'DVDS'  ! variation of sea ice volume due to dynamics in the southern hemisphere
513      titvar(35)= 'DVTN'  ! variation of sea ice volume due to thermo in the   northern hemisphere
514      titvar(36)= 'DVTS'  ! variation of sea ice volume due to thermo in the   southern hemisphere
515      titvar(37)= 'TG1N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 1 
516      titvar(38)= 'TG1S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 1 
517      titvar(39)= 'TG2N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 2 
518      titvar(40)= 'TG2S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 2 
519      titvar(41)= 'TG3N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 3 
520      titvar(42)= 'TG3S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 3 
521      titvar(43)= 'TG4N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 4 
522      titvar(44)= 'TG4S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 4 
523      titvar(45)= 'TG5N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 5 
524      titvar(46)= 'TG5S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 5 
525      titvar(47)= 'LA_N'  ! lateral accretion growth rate, northern hemisphere
526      titvar(48)= 'LA_S'  ! lateral accretion growth rate, southern hemisphere
527      titvar(49)= 'SF_N'  ! Salinity FY, NH
528      titvar(50)= 'SF_S'  ! Salinity FY, SH
529      titvar(51)= 'SF_N'  ! Salinity MY, NH
530      titvar(52)= 'SF_S'  ! Salinity MY, SH
531      titvar(53)= 'Fs_N'  ! Total salt flux NH
532      titvar(54)= 'Fs_S'  ! Total salt flux SH
533      titvar(55)= 'FsbN'  ! Salt - brine drainage flux NH
534      titvar(56)= 'FsbS'  ! Salt - brine drainage flux SH
535      titvar(57)= 'FseN'  ! Salt - Equivalent salt flux NH
536      titvar(58)= 'FseS'  ! Salt - Equivalent salt flux SH
537      titvar(59)= 'SSTN'  ! SST, NH
538      titvar(60)= 'SSTS'  ! SST, SH
539      titvar(61)= 'SSSN'  ! SSS, NH
540      titvar(62)= 'SSSS'  ! SSS, SH
541      titvar(63)= 'TsuN'  ! Tsu, NH
542      titvar(64)= 'TsuS'  ! Tsu, SH
543      titvar(65)= 'TsnN'  ! Tsn, NH
544      titvar(66)= 'TsnS'  ! Tsn, SH
545      titvar(67)= 'ei_N'  ! ei, NH
546      titvar(68)= 'ei_S'  ! ei, SH
547      titvar(69)= 'vi1N'  ! vi1, NH
548      titvar(70)= 'vi1S'  ! vi1, SH
549      titvar(71)= 'vi2N'  ! vi2, NH
550      titvar(72)= 'vi2S'  ! vi2, SH
551      titvar(73)= 'vi3N'  ! vi3, NH
552      titvar(74)= 'vi3S'  ! vi3, SH
553      titvar(75)= 'vi4N'  ! vi4, NH
554      titvar(76)= 'vi4S'  ! vi4, SH
555      titvar(77)= 'vi5N'  ! vi5, NH
556      titvar(78)= 'vi5S'  ! vi5, SH
557      titvar(79)= 'vi6N'  ! vi6, NH
558      titvar(80)= 'vi6S'  ! vi6, SH
559      titvar(81)= 'fmaN'  ! mass flux in the ocean, NH
560      titvar(82)= 'fmaS'  ! mass flux in the ocean, SH
561      titvar(83)= 'AFSE'  ! Fram Strait Area export
562      titvar(84)= 'VFSE'  ! Fram Strait Volume export
563      nvinfo = 84
564
565      ! Definition et Ecriture de l'entete : nombre d'enregistrements
566      ndeb   = ( nstart - 1 ) / ninfo
567      IF( nstart == 1 ) ndeb = -1
568
569      nferme = ( nstart - 1 + nitrun) / ninfo
570      ntot   = nferme - ndeb
571      ndeb   = ninfo * ( 1 + ndeb )
572      nferme = ninfo * nferme
573
574      ! definition of formats
575      WRITE( fmtw  , '(A,I3,A2,I1,A)' )  '(', nfrinf, '(A', jpchsep, ','//fmtinf//'))'
576      WRITE( fmtr  , '(A,I3,A,I1,A)'  )  '(', nfrinf, '(', jpchsep, 'X,'//fmtinf//'))'
577      WRITE( fmtitr, '(A,I3,A,I1,A)'  )  '(', nvinfo, 'A', jpchinf, ')'
578
579      ! opening  "ice_evolu" file
580      IF ( lk_mpp ) THEN
581         WRITE(clname,FMT="('ice.evolu_',I4.4)") narea-1
582      ELSE
583         clname = 'ice.evolu'
584      END IF
585      irecl = ( jpchinf + 1 ) * nvinfo 
586      CALL ctlopn( numevo_ice, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',    &
587         &         irecl, numout, narea, 1 )
588
589      !- ecriture de 2 lignes d''entete :
590      WRITE(numevo_ice,1000) fmtr, fmtw, fmtitr, nvinfo, ntot, 0, nfrinf
591      zxx0 = 0.001 * REAL(ninfo)
592      zxx1 = 0.001 * REAL(ndeb)
593      WRITE(numevo_ice,1111) REAL(jpchinf), 0., zxx1, zxx0, 0., 0., 0
594
595      !- ecriture de 2 lignes de titre :
596      WRITE(numevo_ice,'(A,I8,A,I8,A,I5)')                                      &
597         'Evolution chronologique - Experience '//cexper   &
598         //'   de', ndeb, ' a', nferme, ' pas', ninfo
599      WRITE(numevo_ice,fmtitr) ( titvar(jv), jv = 1, nvinfo )
600
601
602      !--preparation de "titvar" pour l''ecriture parmi les valeurs numeriques :
603      DO  jv = 2 , nvinfo
604         titinf     = titvar(jv)(:jpchinf)
605         titvar(jv) = '  '//titinf
606      END DO
607
608      !--Initialisation of the arrays for the accumulation
609      DO  jv = 1, nvinfo
610         vinfom(jv) = 0.
611      END DO
612      naveg = 0
613
6141000  FORMAT( 3(A20),4(1x,I6) )
6151111  FORMAT( 3(F7.1,1X,F7.3,1X),I3,A ) 
616
617   END SUBROUTINE lim_dia_init
618
619#else
620   !!----------------------------------------------------------------------
621   !!   Default option :                               NO LIM sea-ice model
622   !!----------------------------------------------------------------------
623CONTAINS
624   SUBROUTINE lim_dia         ! Empty routine
625   END SUBROUTINE lim_dia
626#endif
627
628   !!======================================================================
629END MODULE limdia