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.
limdyn.F90 in trunk/NEMO/LIM_SRC – NEMO

source: trunk/NEMO/LIM_SRC/limdyn.F90 @ 106

Last change on this file since 106 was 106, checked in by opalod, 20 years ago

CT : UPDATE067 : Add control indices nictl, njctl used in SUM function output to compare mono versus multi procs runs
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.8 KB
Line 
1MODULE limdyn
2   !!======================================================================
3   !!                     ***  MODULE  limdyn  ***
4   !!   Sea-Ice dynamics : 
5   !!======================================================================
6#if defined key_ice_lim
7   !!----------------------------------------------------------------------
8   !!   'key_ice_lim' :                                   LIM sea-ice model
9   !!----------------------------------------------------------------------
10   !!    lim_dyn      : computes ice velocities
11   !!    lim_dyn_init : initialization and namelist read
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE phycst
15   USE in_out_manager  ! I/O manager
16   USE dom_ice
17   USE dom_oce         ! ocean space and time domain
18   USE ice
19   USE ice_oce
20   USE iceini
21   USE limistate
22   USE limrhg          ! ice rheology
23   USE lbclnk
24   USE lib_mpp
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Accessibility
30   PUBLIC lim_dyn  ! routine called by ice_step
31
32   !! * Module variables
33   REAL(wp)  ::  rone    = 1.e0   ! constant value
34
35   !!----------------------------------------------------------------------
36   !!   LIM 2.0 , UCL-LODYC-IPSL  (2003)
37   !!----------------------------------------------------------------------
38
39CONTAINS
40
41   SUBROUTINE lim_dyn
42      !!-------------------------------------------------------------------
43      !!               ***  ROUTINE lim_dyn  ***
44      !!               
45      !! ** Purpose :   compute ice velocity and ocean-ice stress
46      !!               
47      !! ** Method  :
48      !!
49      !! ** Action  : - Initialisation
50      !!              - Call of the dynamic routine for each hemisphere
51      !!              - computation of the stress at the ocean surface         
52      !!              - treatment of the case if no ice dynamic
53      !! History :
54      !!   1.0  !  01-04  (LIM)  Original code
55      !!   2.0  !  02-08  (C. Ethe, G. Madec)  F90, mpp
56      !!---------------------------------------------------------------------
57      !! * Loal variables
58      INTEGER ::   ji, jj             ! dummy loop indices
59      INTEGER ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology
60      REAL(wp) ::   &
61         ztairx, ztairy,           &  ! tempory scalars
62         zsang , zmod,             &
63         ztglx , ztgly ,           &
64         zt11, zt12, zt21, zt22 ,  &
65         zustm, zsfrld, zsfrldm4,  &
66         zu_ice, zv_ice, ztair2
67      REAL(wp),DIMENSION(jpj) ::   &
68         zind,                     &  ! i-averaged indicator of sea-ice
69         zmsk                         ! i-averaged of tmask
70      !!---------------------------------------------------------------------
71
72      IF( numit == nstart  )   CALL lim_dyn_init   ! Initialization (first time-step only)
73     
74      IF ( ln_limdyn ) THEN
75
76         ! Mean ice and snow thicknesses.         
77         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)
78         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:)
79
80         u_oce(:,:)  = u_io(:,:) * tmu(:,:)
81         v_oce(:,:)  = v_io(:,:) * tmu(:,:)
82       
83         !                                         ! Rheology (ice dynamics)
84         !                                         ! ========
85         
86         !  Define the j-limits where ice rheology is computed
87         ! ---------------------------------------------------
88         
89         IF( lk_mpp ) THEN                    ! mpp: compute over the whole domain
90            i_j1 = 1   
91            i_jpj = jpj
92            IF(l_ctl)   WRITE(numout,*) 'lim_dyn :    i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
93            CALL lim_rhg( i_j1, i_jpj )
94
95         ELSE                                 ! optimization of the computational area
96
97            DO jj = 1, jpj
98               zind(jj) = SUM( frld (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
99               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
100   !!i         write(numout,*) narea, 'limdyn' , jj, zind(jj), zmsk(jj)
101            END DO
102
103            IF( l_jeq ) THEN                     ! local domain include both hemisphere
104               !                                 ! Rheology is computed in each hemisphere
105               !                                 ! only over the ice cover latitude strip
106               ! Northern hemisphere
107               i_j1  = njeq
108               i_jpj = jpj
109               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
110                  i_j1 = i_j1 + 1
111               END DO
112               i_j1 = MAX( 1, i_j1-1 )
113               IF(l_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
114   
115               CALL lim_rhg( i_j1, i_jpj )
116   
117               ! Southern hemisphere
118               i_j1  =  1 
119               i_jpj = njeq
120               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
121                  i_jpj = i_jpj - 1
122               END DO
123               i_jpj = MIN( jpj, i_jpj+2 )
124               IF(l_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
125   
126               CALL lim_rhg( i_j1, i_jpj )
127   
128            ELSE                                 ! local domain extends over one hemisphere only
129               !                                 ! Rheology is computed only over the ice cover
130               !                                 ! latitude strip
131               i_j1  = 1
132               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
133                  i_j1 = i_j1 + 1
134               END DO
135               i_j1 = MAX( 1, i_j1-1 )
136   
137               i_jpj  = jpj
138               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
139                  i_jpj = i_jpj - 1
140               END DO
141               i_jpj = MIN( jpj, i_jpj+2)
142   
143               IF(l_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
144   
145               CALL lim_rhg( i_j1, i_jpj )
146
147            ENDIF
148
149         ENDIF
150
151         u_ice(:,1) = 0.e0       !ibug  est-ce vraiment necessaire?
152         v_ice(:,1) = 0.e0
153
154         IF(l_ctl) THEN
155            WRITE(numout,*) ' lim_dyn  : u_oce ', SUM( u_oce(2:nictl,2:njctl) ), ' v_oce ', SUM( v_oce(2:nictl,2:njctl) )
156            WRITE(numout,*) ' lim_dyn  : u_ice ', SUM( u_ice(2:nictl,2:njctl) ), ' v_ice ', SUM( v_ice(2:nictl,2:njctl) )
157         ENDIF
158         
159         !                                         ! Ice-Ocean stress
160         !                                         ! ================
161         DO jj = 2, jpjm1
162            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
163            DO ji = 2, jpim1
164               ! computation of wind stress over ocean in X and Y direction
165#if defined key_coupled && defined key_lim_cp1
166               ztairx =  frld(ji-1,jj  ) * gtaux(ji-1,jj  ) + frld(ji,jj  ) * gtaux(ji,jj  )      &
167                  &    + frld(ji-1,jj-1) * gtaux(ji-1,jj-1) + frld(ji,jj-1) * gtaux(ji,jj-1)
168
169               ztairy =  frld(ji-1,jj  ) * gtauy(ji-1,jj  ) + frld(ji,jj  ) * gtauy(ji,jj  )      &
170                  &    + frld(ji-1,jj-1) * gtauy(ji-1,jj-1) + frld(ji,jj-1) * gtauy(ji,jj-1)
171#else
172               zsfrld  = frld(ji,jj) + frld(ji-1,jj) + frld(ji-1,jj-1) + frld(ji,jj-1)
173               ztairx  = zsfrld * gtaux(ji,jj)
174               ztairy  = zsfrld * gtauy(ji,jj)
175#endif
176               zsfrldm4 = 4 - frld(ji,jj) - frld(ji-1,jj) - frld(ji-1,jj-1) - frld(ji,jj-1)
177               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
178               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
179               zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 
180               ztglx   = zsfrldm4 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
181               ztgly   = zsfrldm4 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
182
183               tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 4 * rau0 )
184               tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 4 * rau0 )
185            END DO
186         END DO
187         
188         ! computation of friction velocity
189         DO jj = 2, jpjm1
190            DO ji = 2, jpim1
191
192               zu_ice   = u_ice(ji-1,jj-1) - u_oce(ji-1,jj-1)
193               zv_ice   = v_ice(ji-1,jj-1) - v_oce(ji-1,jj-1)
194               zt11  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice )
195
196               zu_ice   = u_ice(ji-1,jj) - u_oce(ji-1,jj)
197               zv_ice   = v_ice(ji-1,jj) - v_oce(ji-1,jj)
198               zt12  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
199
200               zu_ice   = u_ice(ji,jj-1) - u_oce(ji,jj-1)
201               zv_ice   = v_ice(ji,jj-1) - v_oce(ji,jj-1)
202               zt21  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
203
204               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
205               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
206               zt22  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
207
208               ztair2 = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj)
209
210               zustm =  ( 1 - frld(ji,jj) ) * 0.25 * ( zt11 + zt12 + zt21 + zt22 )        &
211                  &  +        frld(ji,jj)   * SQRT( ztair2 )
212
213               ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
214            END DO
215         END DO
216
217       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
218                   
219          DO jj = 2, jpjm1
220             DO ji = 2, jpim1
221#if defined key_coupled && defined key_lim_cp1
222                tio_u(ji,jj) = - (  gtaux(ji  ,jj  ) + gtaux(ji-1,jj  )       &
223                   &              + gtaux(ji-1,jj-1) + gtaux(ji  ,jj-1) ) / ( 4 * rau0 )
224
225                tio_v(ji,jj) = - (  gtauy(ji  ,jj )  + gtauy(ji-1,jj  )       &
226                   &              + gtauy(ji-1,jj-1) + gtauy(ji  ,jj-1) ) / ( 4 * rau0 )
227#else
228                tio_u(ji,jj) = - gtaux(ji,jj) / rau0
229                tio_v(ji,jj) = - gtauy(ji,jj) / rau0 
230#endif
231                ztair2       = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj)
232                zustm        = SQRT( ztair2  )
233
234                ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
235            END DO
236         END DO
237
238      ENDIF
239
240      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
241      CALL lbc_lnk( tio_u, 'I', -1. )   ! I-point (i.e. ice U-V point)
242      CALL lbc_lnk( tio_v, 'I', -1. )   ! I-point (i.e. ice U-V point)
243
244      IF(l_ctl) THEN
245         WRITE(numout,*) ' lim_dyn  : tio_u ', SUM( tio_u(2:nictl,2:njctl) ), ' tio_v ', SUM( tio_v(2:nictl,2:njctl) )
246         WRITE(numout,*) ' lim_dyn  : ust2s ', SUM( ust2s(2:nictl,2:njctl) )
247      ENDIF
248
249   END SUBROUTINE lim_dyn
250
251
252   SUBROUTINE lim_dyn_init
253      !!-------------------------------------------------------------------
254      !!                  ***  ROUTINE lim_dyn_init  ***
255      !!
256      !! ** Purpose : Physical constants and parameters linked to the ice
257      !!      dynamics
258      !!
259      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
260      !!       parameter values called at the first timestep (nit000)
261      !!
262      !! ** input   :   Namelist namicedyn
263      !!
264      !! history :
265      !!  8.5  ! 03-08 (C. Ethe) original code
266      !!-------------------------------------------------------------------
267      NAMELIST/namicedyn/ epsd, alpha,     &
268         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
269         &                c_rhg, etamn, creepl, ecc, ahi0
270      !!-------------------------------------------------------------------
271
272      ! Define the initial parameters
273      ! -------------------------
274
275      ! Read Namelist namicedyn
276      REWIND ( numnam_ice )
277      READ   ( numnam_ice  , namicedyn )
278      IF(lwp) THEN
279         WRITE(numout,*)
280         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
281         WRITE(numout,*) '~~~~~~~~~~~~'
282         WRITE(numout,*) '       tolerance parameter                              epsd   = ', epsd
283         WRITE(numout,*) '       coefficient for semi-implicit coriolis           alpha  = ', alpha
284         WRITE(numout,*) '       diffusion constant for dynamics                  dm     = ', dm
285         WRITE(numout,*) '       number of sub-time steps for relaxation          nbiter = ', nbiter
286         WRITE(numout,*) '       maximum number of iterations for relaxation      nbitdr = ', nbitdr
287         WRITE(numout,*) '       relaxation constant                              om     = ', om
288         WRITE(numout,*) '       maximum value for the residual of relaxation     resl   = ', resl
289         WRITE(numout,*) '       drag coefficient for oceanic stress              cw     = ', cw
290         WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg
291         WRITE(numout,*) '       first bulk-rheology parameter                    pstar  = ', pstar
292         WRITE(numout,*) '       second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
293         WRITE(numout,*) '       minimun value for viscosity                      etamn  = ', etamn
294         WRITE(numout,*) '       creep limit                                      creepl = ', creepl
295         WRITE(numout,*) '       eccentricity of the elliptical yield curve       ecc    = ', ecc
296         WRITE(numout,*) '       horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
297      ENDIF
298
299      usecc2 = 1.0 / ( ecc * ecc )
300      rhoco  = rau0 * cw
301      angvg  = angvg * rad
302      sangvg = SIN( angvg )
303      cangvg = COS( angvg )
304      pstarh = pstar / 2.0
305
306      !  Diffusion coefficients.
307      ahiu(:,:) = ahi0 * umask(:,:,1)
308      ahiv(:,:) = ahi0 * vmask(:,:,1)
309
310   END SUBROUTINE lim_dyn_init
311
312#else
313   !!----------------------------------------------------------------------
314   !!   Default option          Empty module           NO LIM sea-ice model
315   !!----------------------------------------------------------------------
316CONTAINS
317   SUBROUTINE lim_dyn         ! Empty routine
318   END SUBROUTINE lim_dyn
319#endif 
320
321   !!======================================================================
322END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.