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

source: trunk/NEMO/LIM_SRC/icestp.F90 @ 420

Last change on this file since 420 was 420, checked in by opalod, 19 years ago

nemo_v1_update_041 : CT : New functionality: ice damping in buffer zones
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.4 KB
Line 
1MODULE icestp
2   !!======================================================================
3   !!                       ***  MODULE icestp   ***
4   !!   Sea-Ice model : LIM Sea ice model time-stepping
5   !!======================================================================
6#if defined key_ice_lim
7   !!----------------------------------------------------------------------
8   !!   'key_ice_lim' :                                   Lim sea-ice model
9   !!----------------------------------------------------------------------
10   !!   ice_stp       : sea-ice model time-stepping
11   !!----------------------------------------------------------------------
12   USE dom_oce
13   USE oce  ! dynamics and tracers variables
14   USE in_out_manager
15   USE ice_oce         ! ice variables
16   USE flx_oce         ! forcings variables
17   USE dom_ice
18   USE cpl_oce
19   USE daymod
20   USE phycst          ! Define parameters for the routines
21   USE taumod
22   USE ice
23   USE iceini
24   USE ocesbc
25   USE lbclnk
26   USE limdyn
27   USE limtrp
28   USE limthd
29   USE limflx
30   USE limdia
31   USE limwri
32   USE limrst
33   USE limdmp          ! Ice damping
34   USE prtctl          ! Print control
35
36   IMPLICIT NONE
37   PRIVATE
38
39   !! * Routine accessibility
40   PUBLIC ice_stp  ! called by step.F90
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!-----------------------------------------------------
46   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)
47   !! $Header$
48   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
49   !!-----------------------------------------------------
50
51CONTAINS
52
53   SUBROUTINE ice_stp ( kt )
54      !!---------------------------------------------------------------------
55      !!                  ***  ROUTINE ice_stp  ***
56      !!                   
57      !! ** Purpose :   Louvain la Neuve Sea Ice Model time stepping
58      !!
59      !! ** Action  : - call the ice dynamics routine
60      !!              - call the ice advection/diffusion routine
61      !!              - call the ice thermodynamics routine
62      !!              - call the routine that computes mass and
63      !!                heat fluxes at the ice/ocean interface
64      !!              - save the outputs
65      !!              - save the outputs for restart when necessary
66      !!
67      !! History :
68      !!   1.0  !  99-11  (M. Imbard)  Original code
69      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
70      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
71      !!----------------------------------------------------------------------
72      !! * Arguments
73      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
74
75      !! * Local declarations
76      INTEGER   ::   ji, jj   ! dummy loop indices
77
78      REAL(wp) , DIMENSION(jpi,jpj)    :: &
79         zsss_io, zsss2_io, zsss3_io          ! tempory workspaces
80      !!----------------------------------------------------------------------
81
82      IF( kt == nit000 ) THEN
83         IF( lk_cpl ) THEN
84            IF(lwp) WRITE(numout,*)
85            IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)'
86            IF(lwp) WRITE(numout,*) '~~~~~~~   coupled case'
87         ELSE
88            IF(lwp) WRITE(numout,*)
89            IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)' 
90            IF(lwp) WRITE(numout,*) '~~~~~~~   forced case using bulk formulea'
91         ENDIF
92         !  Initialize fluxes fields
93         gtaux(:,:) = 0.e0
94         gtauy(:,:) = 0.e0
95      ENDIF
96
97      ! Temperature , salinity and horizonta wind
98      ! sst_io  and sss_io, u_io and v_io  are initialized at nit000 in limistate.F90 (or limrst.F90) with :
99      !              sst_io = sst_io + (nfice - 1) * (tn(:,:,1)+rt0 )
100      !              sss_io = sss_io + (nfice - 1) * sn(:,:,1)
101      !              u_io  = u_io  + (nfice - 1) * 0.5 * ( un(ji-1,jj  ,1) + un(ji-1,jj-1,1) )
102      !              v_io  = v_io  + (nfice - 1) * 0.5 * ( vn(ji  ,jj-1,1) + vn(ji-1,jj-1,1) )
103      !    cumulate fields
104      !
105      sst_io(:,:) = sst_io(:,:) + tn(:,:,1) + rt0
106      sss_io(:,:) = sss_io(:,:) + sn(:,:,1)
107
108 
109      ! vectors at F-point
110      DO jj = 2, jpj
111         DO ji = fs_2, jpi   ! vector opt.
112            u_io(ji,jj) = u_io(ji,jj) + 0.5 * ( un(ji-1,jj  ,1) + un(ji-1,jj-1,1) )
113            v_io(ji,jj) = v_io(ji,jj) + 0.5 * ( vn(ji  ,jj-1,1) + vn(ji-1,jj-1,1) )
114         END DO
115      END DO
116     
117     
118      IF( MOD( kt-1, nfice ) == 0 ) THEN
119         
120         ! The LIM model is going to be call
121         sst_io(:,:) = sst_io(:,:) / FLOAT( nfice ) * tmask(:,:,1)
122         sss_io(:,:) = sss_io(:,:) / FLOAT( nfice )
123
124         ! stress from ocean U- and V-points to ice U,V point
125         DO jj = 2, jpj
126            DO ji = fs_2, jpi   ! vector opt.
127               gtaux(ji,jj) = 0.5 * ( taux(ji-1,jj  ) + taux(ji-1,jj-1) )
128               gtauy(ji,jj) = 0.5 * ( tauy(ji  ,jj-1) + tauy(ji-1,jj-1) )
129               u_io  (ji,jj) = u_io(ji,jj) / FLOAT( nfice )
130               v_io  (ji,jj) = v_io(ji,jj) / FLOAT( nfice )
131            END DO
132         END DO
133
134         ! lateral boundary condition
135         CALL lbc_lnk( gtaux(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
136         CALL lbc_lnk( gtauy(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
137         CALL lbc_lnk( u_io (:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
138         CALL lbc_lnk( v_io (:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
139
140!!gmbug  in the ocean freezing point computed as :
141!!gm           fzptn (ji,jj) = ( -0.0575 + 1.710523e-3 * SQRT( sn(ji,jj,1) )   &
142!!gm                                     - 2.154996e-4 *       sn(ji,jj,1)   ) * sn(ji,jj,1)   !!   &
143!!gm           !!                        - 7.53e-4 * pressure
144!!gm
145!!gm!bug this is much more accurate and efficient computation
146!!gm       **************************************************
147!!gm freezing point from unesco:
148!!gm     real function tf(s,p)
149!!gm   function to compute the freezing point of seawater
150!!gm
151!!gm   reference: unesco tech. papers in the marine science no. 28. 1978
152!!gm   eighth report jpots
153!!gm   annex 6 freezing point of seawater f.j. millero pp.29-35.
154!!gm
155!!gm  units:
156!!gm         pressure      p          decibars
157!!gm         salinity      s          pss-78
158!!gm         temperature   tf         degrees celsius
159!!gm         freezing pt.
160!!gm************************************************************
161!!gm  checkvalue: tf= -2.588567 deg. c for s=40.0, p=500. decibars
162!!gm     tf=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p
163!!gm     return
164!!gm     end
165!!gm!bug
166
167
168!!gm      DO jj = 1, jpj
169!!gm         DO ji = 1, jpi
170!!gm            tfu(ji,jj)  = (  rt0 + ( - 0.0575                              &
171!!gm               &                     + 1.710523e-3 * SQRT( sss_io(ji,jj) )   &
172!!gm               &                     - 2.154996e-4 *       sss_io(ji,jj)   ) * sss_io(ji,jj)  ) * tms(ji,jj)
173!!gm         END DO
174!!gm      END DO
175!!gm
176      zsss_io (:,:) = SQRT( sss_io(:,:) ) 
177      zsss2_io(:,:) =  sss_io(:,:) *  sss_io(:,:)
178      zsss3_io(:,:) = zsss_io(:,:) * zsss_io(:,:) * zsss_io(:,:)
179
180      DO jj = 1, jpj
181         DO ji = 1, jpi
182            tfu(ji,jj)  = ABS ( rt0 - 0.0575       *   sss_io(ji,jj)   &
183               &                    + 1.710523e-03 * zsss3_io(ji,jj)   &
184               &                    - 2.154996e-04 * zsss2_io(ji,jj) ) * tms(ji,jj)
185         END DO
186      END DO
187   
188         
189
190         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
191            CALL prt_ctl_info('Ice Forcings ')
192            CALL prt_ctl(tab2d_1=qsr_oce ,clinfo1=' qsr_oce  : ', tab2d_2=qsr_ice , clinfo2=' qsr_ice   : ')
193            CALL prt_ctl(tab2d_1=qnsr_oce,clinfo1=' qnsr_oce : ', tab2d_2=qnsr_ice, clinfo2=' qnsr_ice  : ')
194            CALL prt_ctl(tab2d_1=evap    ,clinfo1=' evap     : ')
195            CALL prt_ctl(tab2d_1=tprecip ,clinfo1=' precip   : ', tab2d_2=sprecip , clinfo2=' Snow      : ')
196            CALL prt_ctl(tab2d_1=gtaux   ,clinfo1=' u-stress : ', tab2d_2=gtauy   , clinfo2=' v-stress  : ')
197            CALL prt_ctl(tab2d_1=sst_io  ,clinfo1=' sst      : ', tab2d_2=sss_io  , clinfo2=' sss       : ')
198            CALL prt_ctl(tab2d_1=u_io    ,clinfo1=' u_io     : ', tab2d_2=v_io    , clinfo2=' v_io      : ')
199            CALL prt_ctl(tab2d_1=hsnif   ,clinfo1=' hsnif  1 : ', tab2d_2=hicif   , clinfo2=' hicif     : ')
200            CALL prt_ctl(tab2d_1=frld    ,clinfo1=' frld   1 : ', tab2d_2=sist    , clinfo2=' sist      : ')
201         ENDIF
202
203         ! Ice model call
204         numit = numit + nfice 
205
206         !                                                           !--------------!
207         CALL lim_dyn                                                ! Ice dynamics !   ( rheology/dynamics )
208         !                                                           !--------------!
209         IF(ln_ctl) THEN
210            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  2 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
211            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   2 : ', tab2d_2=sist  , clinfo2=' sist      : ')
212         ENDIF
213
214
215         !                                                           !---------------!
216         CALL lim_trp                                                ! Ice transport !  ( Advection/diffusion )
217         !                                                           !---------------!
218         IF(ln_ctl) THEN
219            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  3 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
220            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   3 : ', tab2d_2=sist  , clinfo2=' sist      : ')
221         ENDIF
222
223         !                                                           !-------------!
224         IF( ln_limdmp ) CALL lim_dmp(kt)                            ! Ice damping !
225         !                                                           !-------------!
226
227         !                                                           !--------------------!
228         CALL lim_thd                                                ! Ice thermodynamics !
229         !                                                           !--------------------!
230         IF(ln_ctl) THEN
231            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  4 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
232            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   4 : ', tab2d_2=sist  , clinfo2=' sist      : ')
233         ENDIF
234
235
236         ! Mass and heat fluxes from ice to ocean
237         !                                                           !------------------------------!
238         CALL lim_flx                                                ! Ice/Ocean Mass & Heat fluxes !
239         !                                                           !------------------------------!
240
241         IF( MOD( numit, ninfo ) == 0 .OR. ntmoy == 1 )  THEN        !-----------------!
242            CALL lim_dia                                             ! Ice Diagnostics !
243         ENDIF                                                       !-----------------!
244
245         !                                                           !-------------!
246         CALL lim_wri                                                ! Ice outputs !
247         !                                                           !-------------!
248
249         IF( MOD( numit, nstock ) == 0 .OR. numit == nlast ) THEN
250            !                                                        !------------------!
251            CALL lim_rst_write( numit )                              ! Ice restart file !
252            !                                                        !------------------!
253         ENDIF
254
255         ! Re-initialization of forcings
256         qsr_oce (:,:) = 0.e0
257         qsr_ice (:,:) = 0.e0
258         qnsr_oce(:,:) = 0.e0
259         qnsr_ice(:,:) = 0.e0 
260         dqns_ice(:,:) = 0.e0 
261         tprecip (:,:) = 0.e0 
262         sprecip (:,:) = 0.e0
263         fr1_i0  (:,:) = 0.e0
264         fr2_i0  (:,:) = 0.e0
265         evap    (:,:) = 0.e0
266#if defined key_coupled 
267         rrunoff (:,:) = 0.e0
268         calving (:,:) = 0.e0
269#else
270         qla_ice (:,:) = 0.e0
271         dqla_ice(:,:) = 0.e0
272#endif
273
274      ENDIF
275
276   END SUBROUTINE ice_stp
277
278#else
279   !!----------------------------------------------------------------------
280   !!   Default option           Dummy module          NO LIM sea-ice model
281   !!----------------------------------------------------------------------
282CONTAINS
283   SUBROUTINE ice_stp ( kt )     ! Dummy routine
284      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
285   END SUBROUTINE ice_stp
286#endif
287
288   !!======================================================================
289END MODULE icestp
Note: See TracBrowser for help on using the repository browser.