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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90 @ 4293

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

small corrections for ice categories hi_max and for shifting ice between categories

  • Property svn:keywords set to Id
File size: 11.5 KB
Line 
1MODULE iceini
2   !!======================================================================
3   !!                       ***  MODULE iceini   ***
4   !!   Sea-ice model : LIM Sea ice model Initialization
5   !!======================================================================
6   !! History :  3.0  ! 2008-03  (M. Vancoppenolle) LIM-3 original code
7   !!            3.3  ! 2010-12  (G. Madec) add call to lim_thd_init and lim_thd_sal_init
8   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                     LIM sea-ice model
13   !!----------------------------------------------------------------------
14   !!   ice_init      : sea-ice model initialization
15   !!----------------------------------------------------------------------
16   USE phycst         ! physical constants
17   USE dom_oce        ! ocean domain
18   USE sbc_oce        ! Surface boundary condition: ocean fields
19   USE sbc_ice        ! Surface boundary condition: ice   fields
20   USE ice            ! LIM variables
21   USE par_ice        ! LIM parameters
22   USE dom_ice        ! LIM domain
23   USE thd_ice        ! LIM thermodynamical variables
24   USE limitd_me      ! LIM ice thickness distribution
25   USE limmsh         ! LIM mesh
26   USE limistate      ! LIM initial state
27   USE limrst         ! LIM restart
28   USE limthd         ! LIM ice thermodynamics
29   USE limthd_sal     ! LIM ice thermodynamics: salinity
30   USE limvar         ! LIM variables
31   USE limsbc         ! LIM surface boundary condition
32   USE in_out_manager ! I/O manager
33   USE lib_mpp        ! MPP library
34   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   ice_init   ! called by sbcice_lim.F90
40
41   !!----------------------------------------------------------------------
42   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE ice_init
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE ice_init  ***
51      !!
52      !! ** purpose :   Allocate all the dynamic arrays of the LIM-3 modules
53      !!----------------------------------------------------------------------
54      INTEGER :: ierr
55      !!----------------------------------------------------------------------
56
57      !                                ! Allocate the ice arrays
58      ierr =        ice_alloc        ()      ! ice variables
59      ierr = ierr + dom_ice_alloc    ()      ! domain
60      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing
61      ierr = ierr + thd_ice_alloc    ()      ! thermodynamics
62      ierr = ierr + lim_itd_me_alloc ()      ! ice thickness distribution - mechanics
63      !
64      IF( lk_mpp    )   CALL mpp_sum( ierr )
65      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
66      !
67      !                                ! adequation jpk versus ice/snow layers/categories
68      IF( jpl   > jpk  .OR.  jpm    > jpk .OR.                                    &
69          jkmax > jpk  .OR.  nlay_s > jpk      )   CALL ctl_stop( 'STOP',         &
70         &     'ice_init: the 3rd dimension of workspace arrays is too small.',   &
71         &     'use more ocean levels or less ice/snow layers/categories.' )
72
73      !                                ! Open the namelist file
74      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
75      !
76      CALL ice_run                     ! set some ice run parameters
77      !
78      CALL lim_thd_init                ! set ice thermodynics parameters
79      !
80      CALL lim_thd_sal_init            ! set ice salinity parameters
81      !
82      rdt_ice   = nn_fsbc * rdttra(1)  ! sea-ice timestep
83      r1_rdtice = 1._wp / rdt_ice      ! sea-ice timestep inverse
84      !
85      CALL lim_msh                     ! ice mesh initialization
86      !
87      CALL lim_itd_ini                 ! ice thickness distribution initialization
88      !
89      !                                ! Initial sea-ice state
90      IF( .NOT.ln_rstart ) THEN              ! start from rest
91         numit = 0
92         numit = nit000 - 1
93         CALL lim_istate                        ! start from rest: sea-ice deduced from sst
94         CALL lim_var_agg(1)                    ! aggregate category variables in bulk variables
95         CALL lim_var_glo2eqv                   ! convert global variables in equivalent variables
96      ELSE                                   ! start from a restart file
97         CALL lim_rst_read                      ! read the restart file
98         numit = nit000 - 1
99         CALL lim_var_agg(1)                    ! aggregate ice variables
100         CALL lim_var_glo2eqv                   ! convert global var in equivalent variables
101      ENDIF
102      !
103      hi_max(jpl) = 99._wp  ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl)
104      !
105      CALL lim_sbc_init                ! ice surface boundary condition   
106      !
107      fr_i(:,:) = at_i(:,:)           ! initialisation of sea-ice fraction
108      tn_ice(:,:,:) = t_su(:,:,:)
109      !
110      nstart = numit  + nn_fsbc     
111      nitrun = nitend - nit000 + 1 
112      nlast  = numit  + nitrun 
113      !
114      IF( nstock == 0  )   nstock = nlast + 1
115      !
116   END SUBROUTINE ice_init
117
118
119   SUBROUTINE ice_run
120      !!-------------------------------------------------------------------
121      !!                  ***  ROUTINE ice_run ***
122      !!                 
123      !! ** Purpose :   Definition some run parameter for ice model
124      !!
125      !! ** Method  :   Read the namicerun namelist and check the parameter
126      !!              values called at the first timestep (nit000)
127      !!
128      !! ** input   :   Namelist namicerun
129      !!-------------------------------------------------------------------
130      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb, ln_limdiaout
131      !!-------------------------------------------------------------------
132      !                   
133      REWIND( numnam_ice )                ! Read Namelist namicerun
134      READ  ( numnam_ice , namicerun )
135      !
136      !IF( lk_mpp .AND. ln_nicep ) THEN
137      !   ln_nicep = .FALSE.
138      !   CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' )
139      !ENDIF
140      !
141      IF(lwp) THEN                        ! control print
142         WRITE(numout,*)
143         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice'
144         WRITE(numout,*) ' ~~~~~~'
145         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn
146         WRITE(numout,*) '   maximum ice concentration                               = ', amax
147         WRITE(numout,*) '   atmospheric drag over sea ice                           = ', cai
148         WRITE(numout,*) '   atmospheric drag over ocean                             = ', cao
149         WRITE(numout,*) '   Several ice points in the ice or not in ocean.output    = ', ln_nicep
150         WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb
151         WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout
152      ENDIF
153      !
154   END SUBROUTINE ice_run
155
156
157   SUBROUTINE lim_itd_ini
158      !!------------------------------------------------------------------
159      !!                ***  ROUTINE lim_itd_ini ***
160      !!
161      !! ** Purpose :   Initializes the ice thickness distribution
162      !! ** Method  :   ...
163      !!    Note    : hi_max(jpl) is here set up to a value close to 7 m for
164      !!              limistate (only) and is changed to 99 m in ice_init
165      !!------------------------------------------------------------------
166      INTEGER  ::   jl, jm               ! dummy loop index
167      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
168      !!------------------------------------------------------------------
169
170      IF(lwp) WRITE(numout,*)
171      IF(lwp) WRITE(numout,*) 'lim_itd_ini : Initialization of ice thickness distribution '
172      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
173
174      !------------------------------------------------------------------------------!
175      ! 1) Ice thickness distribution parameters initialization   
176      !------------------------------------------------------------------------------!
177
178      !- Types boundaries (integer)
179      !----------------------------
180      ice_cat_bounds(1,1) = 1
181      ice_cat_bounds(1,2) = jpl
182
183      !- Number of ice thickness categories in each ice type
184      DO jm = 1, jpm
185         ice_ncat_types(jm) = ice_cat_bounds(jm,2) - ice_cat_bounds(jm,1) + 1 
186      END DO
187
188      !- Make the correspondence between thickness categories and ice types
189      !---------------------------------------------------------------------
190      DO jm = 1, jpm       !over types
191         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) !over thickness categories
192            ice_types(jl) = jm
193         END DO
194      END DO
195
196      IF(lwp) THEN 
197         WRITE(numout,*) ' Number of ice types jpm =      ', jpm
198         WRITE(numout,*) ' Number of ice categories jpl = ', jpl
199         DO jm = 1, jpm
200            WRITE(numout,*) ' Ice type ', jm
201            WRITE(numout,*) ' Number of thickness categories ', ice_ncat_types(jm)
202            WRITE(numout,*) ' Thickness category boundaries  ', ice_cat_bounds(jm,1:2)
203         END DO
204         WRITE(numout,*) 'Ice type vector', ice_types(1:jpl)
205         WRITE(numout,*)
206      ENDIF
207
208      !- Thickness categories boundaries
209      !----------------------------------
210      hi_max(:) = 0._wp
211      hi_max_typ(:,:) = 0._wp
212
213      !- Type 1 - undeformed ice
214      zc1 =  3._wp / REAL( ice_cat_bounds(1,2) - ice_cat_bounds(1,1) + 1 , wp )
215      zc2 = 10._wp * zc1
216      zc3 =  3._wp
217
218      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
219         zx1 = REAL( jl-1 , wp ) / REAL( ice_cat_bounds(1,2) - ice_cat_bounds(1,1) + 1 , wp )
220         hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) )
221      END DO
222
223      !- Fill in the hi_max_typ vector, useful in other circumstances
224      ! Tricky trick: hi_max_typ is actually not used in the code and will be removed in a
225      ! next flyspray at this time, the tricky trick will also be removed (Martin, march 08)
226      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
227         hi_max_typ(jl,1) = hi_max(jl)
228      END DO
229
230      IF(lwp) WRITE(numout,*) ' Thickness category boundaries independently of ice type '
231      IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
232
233      IF(lwp) WRITE(numout,*) ' Thickness category boundaries inside ice types '
234      IF(lwp) THEN
235         DO jm = 1, jpm
236            WRITE(numout,*) ' Type number ', jm
237            WRITE(numout,*) ' hi_max_typ : ', hi_max_typ(0:ice_ncat_types(jm),jm)
238         END DO
239      ENDIF
240      !
241      DO jl = 1, jpl
242         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
243      END DO
244      !
245      !
246   END SUBROUTINE lim_itd_ini
247
248#else
249   !!----------------------------------------------------------------------
250   !!   Default option :        Empty module           NO LIM sea-ice model
251   !!----------------------------------------------------------------------
252CONTAINS
253   SUBROUTINE ice_init        ! Empty routine
254   END SUBROUTINE ice_init
255#endif
256
257   !!======================================================================
258END MODULE iceini
Note: See TracBrowser for help on using the repository browser.