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 in branches/dev_002_LIM/NEMO/LIM_SRC_3 – NEMO

source: branches/dev_002_LIM/NEMO/LIM_SRC_3/limdia.F90 @ 830

Last change on this file since 830 was 825, checked in by ctlod, 16 years ago

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