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 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 10 years ago

remove remaining bugs in LIM3, so that it can run in both regional and global config
  • Property svn:keywords set to Id
File size: 12.4 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 reference and configuration namelist files and namelist output file
74      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
75      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
76      CALL ctl_opn( numoni,         'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
77      !
78      CALL ice_run                     ! set some ice run parameters
79      !
80      CALL lim_thd_init                ! set ice thermodynics parameters
81      !
82      CALL lim_thd_sal_init            ! set ice salinity parameters
83      !
84      rdt_ice   = nn_fsbc * rdttra(1)  ! sea-ice timestep
85      r1_rdtice = 1._wp / rdt_ice      ! sea-ice timestep inverse
86      !
87      CALL lim_msh                     ! ice mesh initialization
88      !
89      CALL lim_itd_ini                 ! ice thickness distribution initialization
90      !
91      !                                ! Initial sea-ice state
92      IF( .NOT. ln_rstart ) THEN              ! start from rest
93         numit = 0
94         numit = nit000 - 1
95         CALL lim_istate                        ! start from rest: sea-ice deduced from sst
96         CALL lim_var_agg(1)                    ! aggregate category variables in bulk variables
97         CALL lim_var_glo2eqv                   ! convert global variables in equivalent variables
98      ELSE                                   ! start from a restart file
99         CALL lim_rst_read                      ! read the restart file
100         numit = nit000 - 1
101         CALL lim_var_agg(1)                    ! aggregate ice variables
102         CALL lim_var_glo2eqv                   ! convert global var in equivalent variables
103      ENDIF
104      !
105      hi_max(jpl) = 99._wp  ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl)
106      !
107      CALL lim_sbc_init                 ! ice surface boundary condition   
108      !
109      fr_i(:,:) = at_i(:,:)             ! initialisation of sea-ice fraction
110      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
111      !
112      nstart = numit  + nn_fsbc     
113      nitrun = nitend - nit000 + 1 
114      nlast  = numit  + nitrun 
115      !
116      IF( nstock == 0  )   nstock = nlast + 1
117      !
118   END SUBROUTINE ice_init
119
120
121   SUBROUTINE ice_run
122      !!-------------------------------------------------------------------
123      !!                  ***  ROUTINE ice_run ***
124      !!                 
125      !! ** Purpose :   Definition some run parameter for ice model
126      !!
127      !! ** Method  :   Read the namicerun namelist and check the parameter
128      !!              values called at the first timestep (nit000)
129      !!
130      !! ** input   :   Namelist namicerun
131      !!-------------------------------------------------------------------
132      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb, ln_limdiaout
133      INTEGER  ::   ios                 ! Local integer output status for namelist read
134      !!-------------------------------------------------------------------
135      !                   
136      REWIND( numnam_ice_ref )              ! Namelist namicerun in reference namelist : Parameters for ice
137      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
138901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
139
140      REWIND( numnam_ice_cfg )              ! Namelist namicerun in configuration namelist : Parameters for ice
141      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
142902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
143      WRITE ( numoni, namicerun )
144      !
145      !IF( lk_mpp .AND. ln_nicep ) THEN
146      !   ln_nicep = .FALSE.
147      !   CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' )
148      !ENDIF
149      !
150      IF(lwp) THEN                        ! control print
151         WRITE(numout,*)
152         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice'
153         WRITE(numout,*) ' ~~~~~~'
154         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn
155         WRITE(numout,*) '   maximum ice concentration                               = ', amax
156         WRITE(numout,*) '   atmospheric drag over sea ice                           = ', cai
157         WRITE(numout,*) '   atmospheric drag over ocean                             = ', cao
158         WRITE(numout,*) '   Several ice points in the ice or not in ocean.output    = ', ln_nicep
159         WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb
160         WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout
161      ENDIF
162      !
163   END SUBROUTINE ice_run
164
165
166   SUBROUTINE lim_itd_ini
167      !!------------------------------------------------------------------
168      !!                ***  ROUTINE lim_itd_ini ***
169      !!
170      !! ** Purpose :   Initializes the ice thickness distribution
171      !! ** Method  :   ...
172      !!    Note    : hi_max(jpl) is here set up to a value close to 7 m for
173      !!              limistate (only) and is changed to 99 m in ice_init
174      !!------------------------------------------------------------------
175      INTEGER  ::   jl, jm               ! dummy loop index
176      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
177      !!------------------------------------------------------------------
178
179      IF(lwp) WRITE(numout,*)
180      IF(lwp) WRITE(numout,*) 'lim_itd_ini : Initialization of ice thickness distribution '
181      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
182
183      !------------------------------------------------------------------------------!
184      ! 1) Ice thickness distribution parameters initialization   
185      !------------------------------------------------------------------------------!
186
187      !- Types boundaries (integer)
188      !----------------------------
189      ice_cat_bounds(1,1) = 1
190      ice_cat_bounds(1,2) = jpl
191
192      !- Number of ice thickness categories in each ice type
193      DO jm = 1, jpm
194         ice_ncat_types(jm) = ice_cat_bounds(jm,2) - ice_cat_bounds(jm,1) + 1 
195      END DO
196
197      !- Make the correspondence between thickness categories and ice types
198      !---------------------------------------------------------------------
199      DO jm = 1, jpm       !over types
200         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) !over thickness categories
201            ice_types(jl) = jm
202         END DO
203      END DO
204
205      IF(lwp) THEN 
206         WRITE(numout,*) ' Number of ice types jpm =      ', jpm
207         WRITE(numout,*) ' Number of ice categories jpl = ', jpl
208         DO jm = 1, jpm
209            WRITE(numout,*) ' Ice type ', jm
210            WRITE(numout,*) ' Number of thickness categories ', ice_ncat_types(jm)
211            WRITE(numout,*) ' Thickness category boundaries  ', ice_cat_bounds(jm,1:2)
212         END DO
213         WRITE(numout,*) 'Ice type vector', ice_types(1:jpl)
214         WRITE(numout,*)
215      ENDIF
216
217      !- Thickness categories boundaries
218      !----------------------------------
219      hi_max(:) = 0._wp
220      hi_max_typ(:,:) = 0._wp
221
222      !- Type 1 - undeformed ice
223      zc1 =  3._wp / REAL( ice_cat_bounds(1,2) - ice_cat_bounds(1,1) + 1 , wp )
224      zc2 = 10._wp * zc1
225      zc3 =  3._wp
226
227      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
228         zx1 = REAL( jl-1 , wp ) / REAL( ice_cat_bounds(1,2) - ice_cat_bounds(1,1) + 1 , wp )
229         hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) )
230      END DO
231
232      !- Fill in the hi_max_typ vector, useful in other circumstances
233      ! Tricky trick: hi_max_typ is actually not used in the code and will be removed in a
234      ! next flyspray at this time, the tricky trick will also be removed (Martin, march 08)
235      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
236         hi_max_typ(jl,1) = hi_max(jl)
237      END DO
238
239      IF(lwp) WRITE(numout,*) ' Thickness category boundaries independently of ice type '
240      IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
241
242      IF(lwp) WRITE(numout,*) ' Thickness category boundaries inside ice types '
243      IF(lwp) THEN
244         DO jm = 1, jpm
245            WRITE(numout,*) ' Type number ', jm
246            WRITE(numout,*) ' hi_max_typ : ', hi_max_typ(0:ice_ncat_types(jm),jm)
247         END DO
248      ENDIF
249      !
250      DO jl = 1, jpl
251         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
252      END DO
253      !
254      !
255   END SUBROUTINE lim_itd_ini
256
257#else
258   !!----------------------------------------------------------------------
259   !!   Default option :        Empty module           NO LIM sea-ice model
260   !!----------------------------------------------------------------------
261CONTAINS
262   SUBROUTINE ice_init        ! Empty routine
263   END SUBROUTINE ice_init
264#endif
265
266   !!======================================================================
267END MODULE iceini
Note: See TracBrowser for help on using the repository browser.