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.
limtrp_2.F90 in branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90 @ 7738

Last change on this file since 7738 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.

  • Property svn:keywords set to Id
File size: 17.7 KB
Line 
1MODULE limtrp_2
2   !!======================================================================
3   !!                       ***  MODULE limtrp_2   ***
4   !! LIM 2.0 transport ice model : sea-ice advection/diffusion
5   !!======================================================================
6   !! History :  LIM  !  2000-01 (UCL)  Original code
7   !!            2.0  !  2001-05 (G. Madec, R. Hordoir) opa norm
8   !!             -   !  2004-01 (G. Madec, C. Ethe)  F90, mpp
9   !!            3.3  !  2009-05  (G. Garric, C. Bricaud) addition of the lim2_evp case
10   !!----------------------------------------------------------------------
11#if defined key_lim2
12   !!----------------------------------------------------------------------
13   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_trp_2      : advection/diffusion process of sea ice
16   !!   lim_trp_init_2 : initialization and namelist read
17   !!----------------------------------------------------------------------
18   USE phycst          ! physical constant
19   USE sbc_oce         ! ocean surface boundary condition
20   USE dom_oce         ! ocean domain
21   USE in_out_manager  ! I/O manager
22   USE dom_ice_2       ! LIM-2 domain
23   USE ice_2           ! LIM-2 variables
24   USE limistate_2     ! LIM-2 initial state
25   USE limadv_2        ! LIM-2 advection
26   USE limhdf_2        ! LIM-2 horizontal diffusion
27   USE lbclnk          ! lateral boundary conditions -- MPP exchanges
28   USE lib_mpp         ! MPP library
29   USE wrk_nemo        ! work arrays
30# if defined key_agrif
31   USE agrif_lim2_interp ! nesting
32# endif
33   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   lim_trp_2   ! called by sbc_ice_lim_2
39
40   REAL(wp), PUBLIC ::   bound      !: boundary condit. (0.0 no-slip, 1.0 free-slip)
41
42   REAL(wp)  ::   epsi06 = 1.e-06   ! constant values
43   REAL(wp)  ::   epsi03 = 1.e-03 
44   REAL(wp)  ::   epsi16 = 1.e-16 
45   REAL(wp)  ::   rzero  = 0.e0   
46   REAL(wp)  ::   rone   = 1.e0
47
48   !! * Substitution
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
52   !! $Id$
53   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55
56CONTAINS
57
58   SUBROUTINE lim_trp_2( kt )
59      !!-------------------------------------------------------------------
60      !!                   ***  ROUTINE lim_trp_2 ***
61      !!                   
62      !! ** purpose : advection/diffusion process of sea ice
63      !!
64      !! ** method  : variables included in the process are scalar,   
65      !!     other values are considered as second order.
66      !!     For advection, a second order Prather scheme is used. 
67      !!
68      !! ** action :
69      !!---------------------------------------------------------------------
70      INTEGER, INTENT(in) ::   kt     ! number of iteration
71      !!
72      INTEGER  ::   ji, jj, jk   ! dummy loop indices
73      INTEGER  ::   initad       ! number of sub-timestep for the advection
74      REAL(wp) ::   zindb  , zindsn , zindic, zacrith   ! local scalars
75      REAL(wp) ::   zusvosn, zusvoic, zignm , zindhe    !   -      -
76      REAL(wp) ::   zvbord , zcfl   , zusnit            !   -      -
77      REAL(wp) ::   zrtt   , ztsn   , ztic1 , ztic2     !   -      -
78      REAL(wp), POINTER, DIMENSION(:,:)  ::   zui_u , zvi_v , zsm             ! 2D workspace
79      REAL(wp), POINTER, DIMENSION(:,:)  ::   zs0ice, zs0sn , zs0a            !  -      -
80      REAL(wp), POINTER, DIMENSION(:,:)  ::   zs0c0 , zs0c1 , zs0c2 , zs0st   !  -      -
81      !---------------------------------------------------------------------
82
83      CALL wrk_alloc( jpi, jpj, zui_u , zvi_v , zsm, zs0ice, zs0sn , zs0a, zs0c0 , zs0c1 , zs0c2 , zs0st )
84
85      IF( kt == nit000  )   CALL lim_trp_init_2      ! Initialization (first time-step only)
86
87# if defined key_agrif
88      CALL agrif_trp_lim2_load      ! First interpolation
89# endif
90
91      zsm(:,:) = area(:,:)
92     
93      IF( ln_limdyn ) THEN
94         !-------------------------------------!
95         !   Advection of sea ice properties   !
96         !-------------------------------------!
97
98         ! ice velocities at ocean U- and V-points (zui_u,zvi_v)
99         ! ---------------------------------------
100         IF( lk_lim2_vp ) THEN      ! VP rheology : B-grid sea-ice dynamics (I-point ice velocity)
101            zvbord = 1._wp + ( 1._wp - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions       
102            DO jj = 1, jpjm1
103               DO ji = 1, jpim1   ! NO vector opt.
104                  zui_u(ji,jj) = ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj)+tmu(ji+1,jj+1), zvbord ) )
105                  zvi_v(ji,jj) = ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji,jj+1)+tmu(ji+1,jj+1), zvbord ) )
106               END DO
107            END DO
108            CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )      ! Lateral boundary conditions
109            !
110         ELSE                       ! EVP rheology : C-grid sea-ice dynamics (u- & v-points ice velocity)
111            zui_u(:,:) = u_ice(:,:)      ! EVP rheology: ice (u,v) at u- and v-points
112            zvi_v(:,:) = v_ice(:,:)
113         ENDIF
114
115         ! CFL test for stability
116         ! ----------------------
117         zcfl  = 0._wp
118         zcfl  = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, :     ) ) * rdt_ice / e1u(1:jpim1, :     ) ) )
119         zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) )
120         !
121         IF(lk_mpp)   CALL mpp_max( zcfl )
122         !
123         IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ', zcfl
124
125         ! content of properties
126         ! ---------------------
127         zs0sn (:,:) =  hsnm(:,:)              * area  (:,:)  ! Snow volume.
128         zs0ice(:,:) =  hicm(:,:)              * area  (:,:)  ! Ice volume.
129         zs0a  (:,:) =  ( 1.0 - frld(:,:) )    * area  (:,:)  ! Surface covered by ice.
130         zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn (:,:)  ! Heat content of the snow layer.
131         zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)  ! Heat content of the first ice layer.
132         zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)  ! Heat content of the second ice layer.
133         zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a  (:,:)  ! Heat reservoir for brine pockets.
134         
135 
136         ! Advection (Prather scheme)
137         ! ---------
138         initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )   ! If ice drift field is too fast,         
139         zusnit = 1.0 / REAL( initad )                              ! split the ice time step in two
140         !
141         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN        !==  odd ice time step:  adv_x then adv_y  ==!
142            DO jk = 1, initad
143               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
144               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
145               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
146               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
147               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
148               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
149               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
150               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
151               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
152               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
153               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
154               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
155               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
156               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
157            END DO
158         ELSE                                                 !==  even ice time step:  adv_x then adv_y  ==!
159            DO jk = 1, initad
160               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
161               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
162               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
163               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
164               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
165               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
166               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
167               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
168               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
169               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
170               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
171               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
172               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
173               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
174            END DO
175         ENDIF
176                       
177         ! recover the properties from their contents
178         ! ------------------------------------------
179!!gm Define in limmsh one for all area = 1 /area  (CPU time saved !)
180         zs0ice(:,:) = zs0ice(:,:) / area(:,:)
181         zs0sn (:,:) = zs0sn (:,:) / area(:,:)
182         zs0a  (:,:) = zs0a  (:,:) / area(:,:)
183         zs0c0 (:,:) = zs0c0 (:,:) / area(:,:)
184         zs0c1 (:,:) = zs0c1 (:,:) / area(:,:)
185         zs0c2 (:,:) = zs0c2 (:,:) / area(:,:)
186         zs0st (:,:) = zs0st (:,:) / area(:,:)
187
188
189         !-------------------------------------!
190         !   Diffusion of sea ice properties   !
191         !-------------------------------------!
192
193         ! Masked eddy diffusivity coefficient at ocean U- and V-points
194         ! ------------------------------------------------------------
195         DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
196            DO ji = 1 , fs_jpim1   ! vector opt.
197               pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji  ,jj) ) ) )   &
198                  &        * ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj) ) ) ) * ahiu(ji,jj)
199               pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji,jj  ) ) ) )   &
200                  &        * ( 1.0 - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1) ) ) ) * ahiv(ji,jj)
201            END DO
202         END DO
203!!gm more readable coding: (and avoid an error in F90 with sign of zero)
204!        DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
205!           DO ji = 1 , fs_jpim1   ! vector opt.
206!              IF( MIN( zs0a(ji,jj) , zs0a(ji+1,jj) ) == 0.e0 )   pahu(ji,jj) = 0._wp
207!              IF( MIN( zs0a(ji,jj) , zs0a(ji,jj+1) ) == 0.e0 )   pahv(ji,jj) = 0._wp
208!           END DO
209!        END DO
210!!gm end
211
212         ! diffusion
213         ! ---------
214         CALL lim_hdf_2( zs0ice )
215         CALL lim_hdf_2( zs0sn  )
216         CALL lim_hdf_2( zs0a   )
217         CALL lim_hdf_2( zs0c0  )
218         CALL lim_hdf_2( zs0c1  )
219         CALL lim_hdf_2( zs0c2  )
220         CALL lim_hdf_2( zs0st  )
221
222!!gm see comment this can be skipped
223         zs0ice(:,:) = MAX( rzero, zs0ice(:,:) * area(:,:) )    !!bug:  useless
224         zs0sn (:,:) = MAX( rzero, zs0sn (:,:) * area(:,:) )    !!bug:  cf /area  just below
225         zs0a  (:,:) = MAX( rzero, zs0a  (:,:) * area(:,:) )    !! caution: the suppression of the 2 changes
226         zs0c0 (:,:) = MAX( rzero, zs0c0 (:,:) * area(:,:) )    !! the last digit of the results
227         zs0c1 (:,:) = MAX( rzero, zs0c1 (:,:) * area(:,:) )
228         zs0c2 (:,:) = MAX( rzero, zs0c2 (:,:) * area(:,:) )
229         zs0st (:,:) = MAX( rzero, zs0st (:,:) * area(:,:) )
230
231
232         !-------------------------------------------------------------------!
233         !   Updating and limitation of sea ice properties after transport   !
234         !-------------------------------------------------------------------!
235         DO jj = 1, jpj
236            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH
237            DO ji = 1, jpi
238               !
239               ! Recover mean values over the grid squares.
240               zs0sn (ji,jj) = MAX( rzero, zs0sn (ji,jj)/area(ji,jj) )
241               zs0ice(ji,jj) = MAX( rzero, zs0ice(ji,jj)/area(ji,jj) )
242               zs0a  (ji,jj) = MAX( rzero, zs0a  (ji,jj)/area(ji,jj) )
243               zs0c0 (ji,jj) = MAX( rzero, zs0c0 (ji,jj)/area(ji,jj) )
244               zs0c1 (ji,jj) = MAX( rzero, zs0c1 (ji,jj)/area(ji,jj) )
245               zs0c2 (ji,jj) = MAX( rzero, zs0c2 (ji,jj)/area(ji,jj) )
246               zs0st (ji,jj) = MAX( rzero, zs0st (ji,jj)/area(ji,jj) )
247
248               ! Recover in situ values.
249               zindb         = MAX( rzero, SIGN( rone, zs0a(ji,jj) - epsi06 ) )
250               zacrith       = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) )
251               zs0a (ji,jj)  = zindb * MIN( zs0a(ji,jj), zacrith )
252               hsnif(ji,jj)  = zindb * ( zs0sn(ji,jj) /MAX( zs0a(ji,jj), epsi16 ) )
253               hicif(ji,jj)  = zindb * ( zs0ice(ji,jj)/MAX( zs0a(ji,jj), epsi16 ) )
254               zindsn        = MAX( rzero, SIGN( rone, hsnif(ji,jj) - epsi06 ) )
255               zindic        = MAX( rzero, SIGN( rone, hicif(ji,jj) - epsi03 ) )
256               zindb         = MAX( zindsn, zindic )
257               zs0a (ji,jj)  = zindb * zs0a(ji,jj)
258               frld (ji,jj)  = 1.0 - zs0a(ji,jj)
259               hsnif(ji,jj)  = zindsn * hsnif(ji,jj)
260               hicif(ji,jj)  = zindic * hicif(ji,jj)
261               zusvosn       = 1.0/MAX( hsnif(ji,jj) * zs0a(ji,jj), epsi16 )
262               zusvoic       = 1.0/MAX( hicif(ji,jj) * zs0a(ji,jj), epsi16 )
263               zignm         = MAX( rzero,  SIGN( rone, hsndif - hsnif(ji,jj) ) )
264               zrtt          = 173.15 * rone 
265               ztsn          =          zignm   * tbif(ji,jj,1)  &
266                              + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) ) 
267               ztic1          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c1(ji,jj) ) , tfu(ji,jj) )
268               ztic2          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c2(ji,jj) ) , tfu(ji,jj) )
269 
270               tbif(ji,jj,1) = zindsn * ztsn  + ( 1.0 - zindsn ) * tfu(ji,jj)               
271               tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj)
272               tbif(ji,jj,3) = zindic * ztic2 + ( 1.0 - zindic ) * tfu(ji,jj)
273               qstoif(ji,jj) = zindb  * xlic * zs0st(ji,jj) /  MAX( zs0a(ji,jj), epsi16 )
274            END DO
275         END DO
276         !
277      ENDIF
278      !
279# if defined key_agrif
280      CALL agrif_trp_lim2      ! Fill boundaries of the fine grid
281# endif
282      !
283      CALL wrk_dealloc( jpi, jpj, zui_u , zvi_v , zsm, zs0ice, zs0sn , zs0a, zs0c0 , zs0c1 , zs0c2 , zs0st )
284      !
285   END SUBROUTINE lim_trp_2
286
287
288   SUBROUTINE lim_trp_init_2
289      !!-------------------------------------------------------------------
290      !!                  ***  ROUTINE lim_trp_init_2  ***
291      !!
292      !! ** Purpose :   initialization of ice advection parameters
293      !!
294      !! ** Method  :   Read the namicetrp namelist and check the parameter
295      !!              values called at the first timestep (nit000)
296      !!
297      !! ** input   :   Namelist namicetrp
298      !!-------------------------------------------------------------------
299      INTEGER  ::   ios                 ! Local integer output status for namelist read
300      NAMELIST/namicetrp/ bound
301      !!-------------------------------------------------------------------
302                   
303      REWIND( numnam_ice_ref )              ! Namelist namicetrp in reference namelist : Ice advection
304      READ  ( numnam_ice_ref, namicetrp, IOSTAT = ios, ERR = 901)
305901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicetrp in reference namelist', lwp )
306
307      REWIND( numnam_ice_cfg )              ! Namelist namicetrp in configuration namelist : Ice advection
308      READ  ( numnam_ice_cfg, namicetrp, IOSTAT = ios, ERR = 902 )
309902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicetrp in configuration namelist', lwp )
310      IF(lwm) WRITE ( numoni, namicetrp )
311
312      IF(lwp) THEN
313         WRITE(numout,*)
314         WRITE(numout,*) 'lim_trp_init_2 : Ice parameters for advection '
315         WRITE(numout,*) '~~~~~~~~~~~~~~'
316         WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound
317      ENDIF
318      !
319   END SUBROUTINE lim_trp_init_2
320
321#else
322   !!----------------------------------------------------------------------
323   !!   Default option         Empty Module                No sea-ice model
324   !!----------------------------------------------------------------------
325CONTAINS
326   SUBROUTINE lim_trp_2        ! Empty routine
327   END SUBROUTINE lim_trp_2
328#endif
329
330   !!======================================================================
331END MODULE limtrp_2
Note: See TracBrowser for help on using the repository browser.