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/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90 @ 5592

Last change on this file since 5592 was 3938, checked in by flavoni, 11 years ago

dev_r3406_CNRS_LIM3: update LIM3, see ticket #1116

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