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

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

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.0 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
25   IMPLICIT NONE
26   PRIVATE
27
28   !! * Accessibility
29   PUBLIC lim_dyn  ! routine called by ice_step
30
31   !! * Module variables
32   REAL(wp)  ::  rone    = 1.0   ! constant value
33
34   !!----------------------------------------------------------------------
35   !!   LIM 2.0 , UCL-LODYC-IPSL  (2003)
36   !!----------------------------------------------------------------------
37
38CONTAINS
39
40   SUBROUTINE lim_dyn
41      !!-------------------------------------------------------------------
42      !!               ***  ROUTINE lim_dyn  ***
43      !!               
44      !! ** Purpose :   compute ice velocity and ocean-ice stress
45      !!               
46      !! ** Method  :
47      !!
48      !! ** Action  : - Initialisation
49      !!              - Call of the dynamic routine for each hemisphere
50      !!              - computation of the stress at the ocean surface         
51      !!              - treatment of the case if no ice dynamic
52      !! History :
53      !!   1.0  !  01-04  (LIM)  Original code
54      !!   2.0  !  02-08  (C. Ethe, G. Madec)  F90, mpp
55      !!---------------------------------------------------------------------
56      !! * Loal variables
57      INTEGER ::   ji, jj,        &  ! dummy loop indices
58         jhemis                      ! jhemis = 1 (NH) ; jhemis = -1 (SH)
59      REAL(wp) ::   &
60         ztairx, ztairy,          &  ! tempory scalars
61         zsang , zmod,            &
62         ztglx , ztgly ,          &
63         zt11, zt12, zt21, zt22 , &
64         zustm, zsfrld, zsfrldm4, &
65         zu_ice, zv_ice, ztair2
66
67      !!---------------------------------------------------------------------
68
69
70      IF( numit == nstart  )   CALL lim_dyn_init   ! Initialization (first time-step only)
71     
72      IF ( ldyn ) THEN
73
74         ! Mean ice and snow thicknesses.         
75         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)
76         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:)
77
78         u_oce(:,:)  = u_io(:,:) * tmu(:,:)
79         v_oce(:,:)  = v_io(:,:) * tmu(:,:)
80       
81         !                                         ! Rheology (ice dynamics)
82         !--  Northern hemisphere                  ! ========
83         jhemis = +1
84         CALL lim_rhg( jhemis )
85
86         !--  Southern hemisphere.
87         jhemis = -1
88         CALL lim_rhg( jhemis )
89
90         u_ice(:,1) = 0.0       !ibug  est-ce vraiment necessaire?
91         v_ice(:,1) = 0.0
92
93         IF( l_ctl .AND. lwp ) THEN
94            WRITE(numout,*) ' lim_dyn  : u_oce ', SUM( u_oce ), ' v_oce ', SUM( v_oce )
95            WRITE(numout,*) ' lim_dyn  : u_ice ', SUM( u_ice ), ' v_ice ', SUM( v_ice )
96         ENDIF
97         
98         !                                         ! Ice-Ocean stress
99         !                                         ! ================
100         DO jj = 2, jpjm1
101            jhemis = SIGN(1, jj - jeq )
102            zsang  = REAL(jhemis) * sangvg
103            DO ji = 2, jpim1
104               ! computation of wind stress over ocean in X and Y direction
105#if defined key_coupled && defined key_lim_cp1
106               ztairx =  frld(ji-1,jj  ) * gtaux(ji-1,jj  ) + frld(ji,jj  ) * gtaux(ji,jj  )      &
107                  &    + frld(ji-1,jj-1) * gtaux(ji-1,jj-1) + frld(ji,jj-1) * gtaux(ji,jj-1)
108
109               ztairy =  frld(ji-1,jj  ) * gtauy(ji-1,jj  ) + frld(ji,jj  ) * gtauy(ji,jj  )      &
110                  &    + frld(ji-1,jj-1) * gtauy(ji-1,jj-1) + frld(ji,jj-1) * gtauy(ji,jj-1)
111#else
112               zsfrld  = frld(ji,jj) + frld(ji-1,jj) + frld(ji-1,jj-1) + frld(ji,jj-1)
113               ztairx  = zsfrld * gtaux(ji,jj)
114               ztairy  = zsfrld * gtauy(ji,jj)
115#endif
116               zsfrldm4 = 4 - frld(ji,jj) - frld(ji-1,jj) - frld(ji-1,jj-1) - frld(ji,jj-1)
117               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
118               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
119               zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 
120               ztglx   = zsfrldm4 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
121               ztgly   = zsfrldm4 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
122
123               tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 4 * rau0 )
124               tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 4 * rau0 )
125            END DO
126         END DO
127         
128         ! computation of friction velocity
129         DO jj = 2, jpjm1
130            DO ji = 2, jpim1
131
132               zu_ice   = u_ice(ji-1,jj-1) - u_oce(ji-1,jj-1)
133               zv_ice   = v_ice(ji-1,jj-1) - v_oce(ji-1,jj-1)
134               zt11  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice )
135
136               zu_ice   = u_ice(ji-1,jj) - u_oce(ji-1,jj)
137               zv_ice   = v_ice(ji-1,jj) - v_oce(ji-1,jj)
138               zt12  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
139
140               zu_ice   = u_ice(ji,jj-1) - u_oce(ji,jj-1)
141               zv_ice   = v_ice(ji,jj-1) - v_oce(ji,jj-1)
142               zt21  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
143
144               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
145               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
146               zt22  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
147
148               ztair2 = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj)
149
150               zustm =  ( 1 - frld(ji,jj) ) * 0.25 * ( zt11 + zt12 + zt21 + zt22 )        &
151                  &  +        frld(ji,jj)   * SQRT( ztair2 )
152
153               ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
154            END DO
155         END DO
156
157       ELSE      ! If no ice dynamics 
158                   
159          DO jj = 2, jpjm1
160             DO ji = 2, jpim1
161#if defined key_coupled && defined key_lim_cp1
162                tio_u(ji,jj) = - (  gtaux(ji  ,jj  ) + gtaux(ji-1,jj  )       &
163                   &              + gtaux(ji-1,jj-1) + gtaux(ji  ,jj-1) ) / ( 4 * rau0 )
164
165                tio_v(ji,jj) = - (  gtauy(ji  ,jj )  + gtauy(ji-1,jj  )       &
166                   &              + gtauy(ji-1,jj-1) + gtauy(ji  ,jj-1) ) / ( 4 * rau0 )
167#else
168                tio_u(ji,jj) = - gtaux(ji,jj) / rau0
169                tio_v(ji,jj) = - gtauy(ji,jj) / rau0 
170#endif
171                ztair2       = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj)
172                zustm        = SQRT( ztair2  )
173
174                ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
175            END DO
176         END DO
177
178      ENDIF
179
180      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
181      CALL lbc_lnk( tio_u, 'I', -1. )   ! I-point (i.e. ice U-V point)
182      CALL lbc_lnk( tio_v, 'I', -1. )   ! I-point (i.e. ice U-V point)
183
184      IF( l_ctl .AND. lwp ) THEN
185         WRITE(numout,*) ' lim_dyn  : tio_u ', SUM( tio_u ), ' tio_v ', SUM( tio_v )
186         WRITE(numout,*) ' lim_dyn  : ust2s ', SUM( ust2s )
187      ENDIF
188
189   END SUBROUTINE lim_dyn
190
191    SUBROUTINE lim_dyn_init
192      !!-------------------------------------------------------------------
193      !!                  ***  ROUTINE lim_dyn_init  ***
194      !!
195      !! ** Purpose : Physical constants and parameters linked to the ice
196      !!      dynamics
197      !!
198      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
199      !!       parameter values called at the first timestep (nit000)
200      !!
201      !! ** input   :   Namelist namicedyn
202      !!
203      !! history :
204      !!  8.5  ! 03-08 (C. Ethe) original code
205      !!-------------------------------------------------------------------
206      NAMELIST/namicedyn/ epsd, nlminn, nlmaxs, alpha,     &
207         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
208         &                c_rhg, etamn, creepl, ecc, ahi0
209      !!-------------------------------------------------------------------
210
211      ! Define the initial parameters
212      ! -------------------------
213
214      ! Read Namelist namicedyn
215      REWIND ( numnam_ice )
216      READ   ( numnam_ice  , namicedyn )
217      IF(lwp) THEN
218         WRITE(numout,*)
219         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
220         WRITE(numout,*) '~~~~~~~~~~~~'
221         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
222         WRITE(numout,*) '   northern minimum index for ice drift             nlminn = ', nlminn
223         WRITE(numout,*) '   southern minimum index for ice drift             nlmaxs = ', nlmaxs
224         WRITE(numout,*) '   coefficient for semi-implicit coriolis           alpha  = ', alpha
225         WRITE(numout,*) '   diffusion constant for dynamics                  dm     = ', dm
226         WRITE(numout,*) '   number of sub-time steps for relaxation          nbiter = ', nbiter
227         WRITE(numout,*) '   maximum number of iterations for relaxation      nbitdr = ', nbitdr
228         WRITE(numout,*) '   relaxation constant                              om     = ', om
229         WRITE(numout,*) '   maximum value for the residual of relaxation     resl   = ', resl
230         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
231         WRITE(numout,*) '   turning angle for oceanic stress                 angvg  = ', angvg
232         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
233         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
234         WRITE(numout,*) '   minimun value for viscosity                      etamn  = ', etamn
235         WRITE(numout,*) '   creep limit                                      creepl = ', creepl
236         WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
237         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
238      ENDIF
239
240      usecc2 = 1.0 / ( ecc * ecc )
241      rhoco  = rau0 * cw
242      angvg  = angvg * rad
243      sangvg = SIN( angvg )
244      cangvg = COS( angvg )
245      pstarh = pstar / 2.0
246
247      !  Diffusion coefficients.
248      ahiu(:,:) = ahi0 * umask(:,:,1)
249      ahiv(:,:) = ahi0 * vmask(:,:,1)
250
251   END SUBROUTINE lim_dyn_init
252
253#else
254   !!----------------------------------------------------------------------
255   !!   Default option          Empty module           NO LIM sea-ice model
256   !!----------------------------------------------------------------------
257CONTAINS
258   SUBROUTINE lim_dyn         ! Empty routine
259   END SUBROUTINE lim_dyn
260#endif 
261
262   !!======================================================================
263END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.