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.
limistate.F90 in branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 4045

Last change on this file since 4045 was 4045, checked in by clem, 9 years ago

new LIM3

  • Property svn:keywords set to Id
File size: 25.6 KB
Line 
1MODULE limistate
2   !!======================================================================
3   !!                     ***  MODULE  limistate  ***
4   !!              Initialisation of diagnostics ice variables
5   !!======================================================================
6   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
7   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
8   !!             -   ! 2012    (C. Rousset) add par_oce (for jp_sal)...bug?
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3' :                                    LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_istate      :  Initialisation of diagnostics ice variables
15   !!   lim_istate_init :  initialization of ice state and namelist read
16   !!----------------------------------------------------------------------
17   USE phycst           ! physical constant
18   USE oce              ! dynamics and tracers variables
19   USE dom_oce          ! ocean domain
20   USE sbc_oce          ! Surface boundary condition: ocean fields
21   USE eosbn2           ! equation of state
22   USE ice              ! sea-ice variables
23   USE par_ice          ! ice parameters
24   USE par_oce          ! ocean parameters
25   USE dom_ice          ! sea-ice domain
26   USE in_out_manager   ! I/O manager
27   USE lbclnk           ! lateral boundary condition - MPP exchanges
28   USE lib_mpp          ! MPP library
29   USE wrk_nemo         ! work arrays
30   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   lim_istate      ! routine called by lim_init.F90
36
37   !                                  !!** init namelist (namiceini) **
38   REAL(wp) ::   ttest    = 2.0_wp     ! threshold water temperature for initial sea ice
39   REAL(wp) ::   hninn    = 0.5_wp     ! initial snow thickness in the north
40   REAL(wp) ::   hginn_u  = 2.5_wp     ! initial ice thickness in the north
41   REAL(wp) ::   aginn_u  = 0.7_wp     ! initial leads area in the north
42   REAL(wp) ::   hginn_d  = 5.0_wp     ! initial ice thickness in the north
43   REAL(wp) ::   aginn_d  = 0.25_wp    ! initial leads area in the north
44   REAL(wp) ::   hnins    = 0.1_wp     ! initial snow thickness in the south
45   REAL(wp) ::   hgins_u  = 1.0_wp     ! initial ice thickness in the south
46   REAL(wp) ::   agins_u  = 0.7_wp     ! initial leads area in the south
47   REAL(wp) ::   hgins_d  = 2.0_wp     ! initial ice thickness in the south
48   REAL(wp) ::   agins_d  = 0.2_wp     ! initial leads area in the south
49   REAL(wp) ::   sinn     = 6.301_wp   ! initial salinity
50   REAL(wp) ::   sins     = 6.301_wp   !
51
52   !!----------------------------------------------------------------------
53   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)
54   !! $Id$
55   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE lim_istate
60      !!-------------------------------------------------------------------
61      !!                    ***  ROUTINE lim_istate  ***
62      !!
63      !! ** Purpose :   defined the sea-ice initial state
64      !!
65      !! ** Method  :   restart from a state defined in a binary file
66      !!                or from arbitrary sea-ice conditions
67      !!-------------------------------------------------------------------
68      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
69      REAL(wp) ::   zeps6, zeps, ztmelts, epsi06   ! local scalars
70      REAL(wp) ::   zvol, zare, zh, zh1, zh2, zh3, zan, zbn, zas, zbs 
71      REAL(wp), POINTER, DIMENSION(:)   ::   zgfactorn, zhin 
72      REAL(wp), POINTER, DIMENSION(:)   ::   zgfactors, zhis
73      REAL(wp), POINTER, DIMENSION(:,:) ::   zidto      ! ice indicator
74      !--------------------------------------------------------------------
75
76      CALL wrk_alloc( jpm, zgfactorn, zgfactors, zhin, zhis )
77      CALL wrk_alloc( jpi, jpj, zidto )
78
79      !--------------------------------------------------------------------
80      ! 1) Preliminary things
81      !--------------------------------------------------------------------
82      epsi06 = 1.e-6_wp
83
84      CALL lim_istate_init     !  reading the initials parameters of the ice
85
86!!gm  in lim2  the initialisation if only done if required in the namelist :
87!!gm      IF( .NOT. ln_limini ) THEN
88!!gm  this should be added in lim3 namelist...
89
90      !--------------------------------------------------------------------
91      ! 2) Ice initialization (hi,hs,frld,t_su,sm_i,t_i,t_s)              |
92      !--------------------------------------------------------------------
93
94      IF(lwp) WRITE(numout,*)
95      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
96      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
97
98      t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius]
99
100      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
101         DO ji = 1, jpi
102            IF( tsn(ji,jj,1,jp_tem)  - t_bo(ji,jj) >= ttest ) THEN   ;   zidto(ji,jj) = 0.e0      ! no ice
103            ELSE                                                     ;   zidto(ji,jj) = 1.e0      !    ice
104            ENDIF
105         END DO
106      END DO
107
108      t_bo(:,:) = t_bo(:,:) + rt0                          ! t_bo converted from Celsius to Kelvin (rt0 over land)
109
110      ! constants for heat contents
111      zeps   = 1.e-20_wp
112      zeps6  = 1.e-06_wp
113
114      ! zgfactor for initial ice distribution
115      zgfactorn(:) = 0._wp
116      zgfactors(:) = 0._wp
117
118      ! first ice type
119      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
120         zhin (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
121         zgfactorn(1) = zgfactorn(1) + exp(-(zhin(1)-hginn_u)*(zhin(1)-hginn_u) * 0.5_wp )
122         zhis (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
123         zgfactors(1) = zgfactors(1) + exp(-(zhis(1)-hgins_u)*(zhis(1)-hgins_u) * 0.5_wp )
124      END DO ! jl
125      zgfactorn(1) = aginn_u / zgfactorn(1)
126      zgfactors(1) = agins_u / zgfactors(1)
127
128      ! -------------
129      ! new distribution, polynom of second order, conserving area and volume
130      zh1 = 0._wp
131      zh2 = 0._wp
132      zh3 = 0._wp
133      DO jl = 1, jpl
134         zh = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
135         zh1 = zh1 + zh
136         zh2 = zh2 + zh * zh
137         zh3 = zh3 + zh * zh * zh
138      END DO
139      IF(lwp) WRITE(numout,*) ' zh1 : ', zh1
140      IF(lwp) WRITE(numout,*) ' zh2 : ', zh2
141      IF(lwp) WRITE(numout,*) ' zh3 : ', zh3
142
143      zvol = aginn_u * hginn_u
144      zare = aginn_u
145      IF( jpl >= 2 ) THEN
146         zbn = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3)
147         zan = ( zare - zbn*zh1 ) / zh2
148      ENDIF
149
150      IF(lwp) WRITE(numout,*) ' zvol: ', zvol
151      IF(lwp) WRITE(numout,*) ' zare: ', zare
152      IF(lwp) WRITE(numout,*) ' zbn : ', zbn 
153      IF(lwp) WRITE(numout,*) ' zan : ', zan 
154
155      zvol = agins_u * hgins_u
156      zare = agins_u
157      IF( jpl >= 2 ) THEN
158         zbs = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3)
159         zas = ( zare - zbs*zh1 ) / zh2
160      ENDIF
161
162      IF(lwp) WRITE(numout,*) ' zvol: ', zvol
163      IF(lwp) WRITE(numout,*) ' zare: ', zare
164      IF(lwp) WRITE(numout,*) ' zbn : ', zbn 
165      IF(lwp) WRITE(numout,*) ' zan : ', zan 
166
167      !end of new lines
168      ! -------------
169!!!
170      ! retour a LIMA_MEC
171      !     ! second ice type
172      !     zdummy  = hi_max(ice_cat_bounds(2,1)-1)
173      !     hi_max(ice_cat_bounds(2,1)-1) = 0.0
174
175      !     ! here to change !!!!
176      !     jm = 2
177      !     DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
178      !        zhin (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
179      !        zhin (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + &
180      !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0
181      !        zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0)
182      !        zhis (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
183      !        zhis (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + &
184      !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0
185      !        zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0)
186      !     END DO ! jl
187      !     zgfactorn(2) = aginn_d / zgfactorn(2)
188      !     zgfactors(2) = agins_d / zgfactors(2)
189      !     hi_max(ice_cat_bounds(2,1)-1) = zdummy
190      ! END retour a LIMA_MEC
191!!!
192
193!!gm  optimisation :  loop over the ice categories inside the ji, jj loop !!!
194
195      DO jj = 1, jpj
196         DO ji = 1, jpi
197
198            !--- Northern hemisphere
199            !----------------------------------------------------------------
200            IF( fcor(ji,jj) >= 0._wp ) THEN   
201
202               !-----------------------
203               ! Ice area / thickness
204               !-----------------------
205
206               IF ( jpl .EQ. 1) THEN ! one category
207
208                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
209                     a_i(ji,jj,jl)    = zidto(ji,jj) * aginn_u
210                     ht_i(ji,jj,jl)   = zidto(ji,jj) * hginn_u
211                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
212                  END DO
213
214               ELSE ! several categories
215
216                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
217                     zhin(1)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
218                     a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(1) * exp(-(zhin(1)-hginn_u)* & 
219                        (zhin(1)-hginn_u)/2.0) , epsi06)
220                     ! new line
221                     a_i(ji,jj,jl)    = zidto(ji,jj) * ( zan * zhin(1) * zhin(1) + zbn * zhin(1) )
222                     ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(1) 
223                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
224                  END DO
225
226               ENDIF
227
228
229!!!
230               ! retour a LIMA_MEC
231               !              !ridged ice
232               !              zdummy  = hi_max(ice_cat_bounds(2,1)-1)
233               !              hi_max(ice_cat_bounds(2,1)-1) = 0.0
234               !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories
235               !                 zhin(2)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
236               !                 a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* &
237               !                                    (zhin(2)-hginn_d)/2.0) , epsi06)
238               !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(2)
239               !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
240               !              END DO
241               !              hi_max(ice_cat_bounds(2,1)-1) = zdummy
242
243               !              !rafted ice
244               !              jl = 6
245               !              a_i(ji,jj,jl)       = 0.0
246               !              ht_i(ji,jj,jl)      = 0.0
247               !              v_i(ji,jj,jl)       = 0.0
248               ! END retour a LIMA_MEC
249!!!
250
251               DO jl = 1, jpl
252
253                  !-------------
254                  ! Snow depth
255                  !-------------
256                  ht_s(ji,jj,jl)   = zidto(ji,jj) * hninn
257                  v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl)
258
259                  !---------------
260                  ! Ice salinity
261                  !---------------
262                  sm_i(ji,jj,jl)   = zidto(ji,jj) * sinn  + ( 1.0 - zidto(ji,jj) ) * 0.1
263                  smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
264
265                  !----------
266                  ! Ice age
267                  !----------
268                  o_i(ji,jj,jl)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) )
269                  oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl)
270
271                  !------------------------------
272                  ! Sea ice surface temperature
273                  !------------------------------
274
275                  t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj)
276
277                  !------------------------------------
278                  ! Snow temperature and heat content
279                  !------------------------------------
280
281                  DO jk = 1, nlay_s
282                     t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt
283                     ! Snow energy of melting
284                     e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
285                     ! Change dimensions
286                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
287                     ! Multiply by volume, so that heat content in 10^9 Joules
288                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
289                        v_s(ji,jj,jl)  / nlay_s
290                  END DO !jk
291
292                  !-----------------------------------------------
293                  ! Ice salinities, temperature and heat content
294                  !-----------------------------------------------
295
296                  DO jk = 1, nlay_i
297                     t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
298                     s_i(ji,jj,jk,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1
299                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
300
301                     ! heat content per unit volume
302                     e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * &
303                        (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
304                        +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &
305                        - rcp      * ( ztmelts - rtt ) &
306                        )
307
308                     ! Correct dimensions to avoid big values
309                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
310
311                     ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
312                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 
313                        area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / &
314                        nlay_i
315                  END DO ! jk
316
317               END DO ! jl
318
319            ELSE ! on fcor
320
321               !--- Southern hemisphere
322               !----------------------------------------------------------------
323
324               !-----------------------
325               ! Ice area / thickness
326               !-----------------------
327
328               IF ( jpl .EQ. 1) THEN ! one category
329
330                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
331                     a_i(ji,jj,jl)    = zidto(ji,jj) * agins_u
332                     ht_i(ji,jj,jl)   = zidto(ji,jj) * hgins_u
333                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
334                  END DO
335
336               ELSE ! several categories
337
338                  !level ice
339                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !over thickness categories
340
341                     zhis(1)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
342                     a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactors(1) * exp(-(zhis(1)-hgins_u) * & 
343                        (zhis(1)-hgins_u)/2.0) , epsi06 )
344                     ! new line square distribution volume conserving
345                     a_i(ji,jj,jl)    = zidto(ji,jj) * ( zas * zhis(1) * zhis(1) + zbs * zhis(1) )
346                     ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(1) 
347                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
348
349                  END DO ! jl
350
351               ENDIF
352
353!!!
354               ! retour a LIMA_MEC
355               !              !ridged ice
356               !              zdummy  = hi_max(ice_cat_bounds(2,1)-1)
357               !              hi_max(ice_cat_bounds(2,1)-1) = 0.0
358               !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories
359               !                 zhis(2)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
360               !                 a_i(ji,jj,jl) = zidto(ji,jj)*MAX( zgfactors(2)   &
361               !                    &          * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 )
362               !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(2)
363               !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
364               !              END DO
365               !              hi_max(ice_cat_bounds(2,1)-1) = zdummy
366
367               !              !rafted ice
368               !              jl = 6
369               !              a_i(ji,jj,jl)       = 0.0
370               !              ht_i(ji,jj,jl)      = 0.0
371               !              v_i(ji,jj,jl)       = 0.0
372               ! END retour a LIMA_MEC
373!!!
374
375               DO jl = 1, jpl !over thickness categories
376
377                  !---------------
378                  ! Snow depth
379                  !---------------
380
381                  ht_s(ji,jj,jl)   = zidto(ji,jj) * hnins
382                  v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl)
383
384                  !---------------
385                  ! Ice salinity
386                  !---------------
387
388                  sm_i(ji,jj,jl)   = zidto(ji,jj) * sins  + ( 1.0 - zidto(ji,jj) ) * 0.1
389                  smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
390
391                  !----------
392                  ! Ice age
393                  !----------
394
395                  o_i(ji,jj,jl)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) )
396                  oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl)
397
398                  !------------------------------
399                  ! Sea ice surface temperature
400                  !------------------------------
401
402                  t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj)
403
404                  !----------------------------------
405                  ! Snow temperature / heat content
406                  !----------------------------------
407
408                  DO jk = 1, nlay_s
409                     t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt
410                     ! Snow energy of melting
411                     e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
412                     ! Change dimensions
413                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
414                     ! Multiply by volume, so that heat content in 10^9 Joules
415                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
416                        v_s(ji,jj,jl)  / nlay_s
417                  END DO
418
419                  !---------------------------------------------
420                  ! Ice temperature, salinity and heat content
421                  !---------------------------------------------
422
423                  DO jk = 1, nlay_i
424                     t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
425                     s_i(ji,jj,jk,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1
426                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
427
428                     ! heat content per unit volume
429                     e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * &
430                        (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
431                        +   lfus  * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &
432                        - rcp      * ( ztmelts - rtt ) &
433                        )
434
435                     ! Correct dimensions to avoid big values
436                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
437
438                     ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
439                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 
440                        area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / &
441                        nlay_i
442                  END DO !jk
443
444               END DO ! jl
445
446            ENDIF ! on fcor
447
448         END DO
449      END DO
450
451      !--------------------------------------------------------------------
452      ! 3) Global ice variables for output diagnostics                    |
453      !--------------------------------------------------------------------
454
455      fsbbq (:,:)     = 0.e0
456      u_ice (:,:)     = 0.e0
457      v_ice (:,:)     = 0.e0
458      stress1_i(:,:)  = 0.0
459      stress2_i(:,:)  = 0.0
460      stress12_i(:,:) = 0.0
461
462      !--------------------------------------------------------------------
463      ! 4) Moments for advection
464      !--------------------------------------------------------------------
465
466      sxopw (:,:) = 0.e0 
467      syopw (:,:) = 0.e0 
468      sxxopw(:,:) = 0.e0 
469      syyopw(:,:) = 0.e0 
470      sxyopw(:,:) = 0.e0
471
472      sxice (:,:,:)  = 0.e0   ;   sxsn (:,:,:)  = 0.e0   ;   sxa  (:,:,:)  = 0.e0
473      syice (:,:,:)  = 0.e0   ;   sysn (:,:,:)  = 0.e0   ;   sya  (:,:,:)  = 0.e0
474      sxxice(:,:,:)  = 0.e0   ;   sxxsn(:,:,:)  = 0.e0   ;   sxxa (:,:,:)  = 0.e0
475      syyice(:,:,:)  = 0.e0   ;   syysn(:,:,:)  = 0.e0   ;   syya (:,:,:)  = 0.e0
476      sxyice(:,:,:)  = 0.e0   ;   sxysn(:,:,:)  = 0.e0   ;   sxya (:,:,:)  = 0.e0
477
478      sxc0  (:,:,:)  = 0.e0   ;   sxe  (:,:,:,:)= 0.e0   
479      syc0  (:,:,:)  = 0.e0   ;   sye  (:,:,:,:)= 0.e0   
480      sxxc0 (:,:,:)  = 0.e0   ;   sxxe (:,:,:,:)= 0.e0   
481      syyc0 (:,:,:)  = 0.e0   ;   syye (:,:,:,:)= 0.e0   
482      sxyc0 (:,:,:)  = 0.e0   ;   sxye (:,:,:,:)= 0.e0   
483
484      sxsal  (:,:,:)  = 0.e0
485      sysal  (:,:,:)  = 0.e0
486      sxxsal (:,:,:)  = 0.e0
487      syysal (:,:,:)  = 0.e0
488      sxysal (:,:,:)  = 0.e0
489
490      !--------------------------------------------------------------------
491      ! 5) Lateral boundary conditions                                    |
492      !--------------------------------------------------------------------
493
494      DO jl = 1, jpl
495         CALL lbc_lnk( a_i(:,:,jl)  , 'T', 1. )
496         CALL lbc_lnk( v_i(:,:,jl)  , 'T', 1. )
497         CALL lbc_lnk( v_s(:,:,jl)  , 'T', 1. )
498         CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. )
499         CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. )
500         !
501         CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. )
502         CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. )
503         CALL lbc_lnk( sm_i(:,:,jl) , 'T', 1. )
504         CALL lbc_lnk( o_i(:,:,jl)  , 'T', 1. )
505         CALL lbc_lnk( t_su(:,:,jl) , 'T', 1. )
506         DO jk = 1, nlay_s
507            CALL lbc_lnk(t_s(:,:,jk,jl), 'T', 1. )
508            CALL lbc_lnk(e_s(:,:,jk,jl), 'T', 1. )
509         END DO
510         DO jk = 1, nlay_i
511            CALL lbc_lnk(t_i(:,:,jk,jl), 'T', 1. )
512            CALL lbc_lnk(e_i(:,:,jk,jl), 'T', 1. )
513         END DO
514         !
515         a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl)
516      END DO
517
518      CALL lbc_lnk( at_i , 'T', 1. )
519      at_i(:,:) = tms(:,:) * at_i(:,:)                       ! put 0 over land
520      !
521      CALL lbc_lnk( fsbbq  , 'T', 1. )
522      !
523      CALL wrk_dealloc( jpm, zgfactorn, zgfactors, zhin, zhis )
524      CALL wrk_dealloc( jpi, jpj, zidto )
525      !
526   END SUBROUTINE lim_istate
527
528
529   SUBROUTINE lim_istate_init
530      !!-------------------------------------------------------------------
531      !!                   ***  ROUTINE lim_istate_init  ***
532      !!       
533      !! ** Purpose : Definition of initial state of the ice
534      !!
535      !! ** Method :   Read the namiceini namelist and check the parameter
536      !!             values called at the first timestep (nit000)
537      !!
538      !! ** input  :   namelist namiceini
539      !!-----------------------------------------------------------------------------
540      NAMELIST/namiceini/ ttest, hninn, hginn_u, aginn_u, hginn_d, aginn_d, hnins,   &
541         &                hgins_u, agins_u, hgins_d, agins_d, sinn, sins
542      !!-----------------------------------------------------------------------------
543      !
544      REWIND ( numnam_ice )               ! Read Namelist namiceini
545      READ   ( numnam_ice , namiceini )
546      !
547      IF(lwp) THEN                        ! control print
548         WRITE(numout,*)
549         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
550         WRITE(numout,*) '~~~~~~~~~~~~~~~'
551         WRITE(numout,*) '   threshold water temp. for initial sea-ice    ttest      = ', ttest
552         WRITE(numout,*) '   initial snow thickness in the north          hninn      = ', hninn
553         WRITE(numout,*) '   initial undef ice thickness in the north     hginn_u    = ', hginn_u
554         WRITE(numout,*) '   initial undef ice concentr. in the north     aginn_u    = ', aginn_u         
555         WRITE(numout,*) '   initial  def  ice thickness in the north     hginn_d    = ', hginn_d
556         WRITE(numout,*) '   initial  def  ice concentr. in the north     aginn_d    = ', aginn_d         
557         WRITE(numout,*) '   initial snow thickness in the south          hnins      = ', hnins 
558         WRITE(numout,*) '   initial undef ice thickness in the north     hgins_u    = ', hgins_u
559         WRITE(numout,*) '   initial undef ice concentr. in the north     agins_u    = ', agins_u         
560         WRITE(numout,*) '   initial  def  ice thickness in the north     hgins_d    = ', hgins_d
561         WRITE(numout,*) '   initial  def  ice concentr. in the north     agins_d    = ', agins_d         
562         WRITE(numout,*) '   initial  ice salinity       in the north     sinn       = ', sinn
563         WRITE(numout,*) '   initial  ice salinity       in the south     sins       = ', sins
564      ENDIF
565      !
566   END SUBROUTINE lim_istate_init
567
568#else
569   !!----------------------------------------------------------------------
570   !!   Default option :         Empty module          NO LIM sea-ice model
571   !!----------------------------------------------------------------------
572CONTAINS
573   SUBROUTINE lim_istate          ! Empty routine
574   END SUBROUTINE lim_istate
575#endif
576
577   !!======================================================================
578END MODULE limistate
Note: See TracBrowser for help on using the repository browser.