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

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

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