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

Last change on this file since 508 was 508, checked in by opalod, 18 years ago

nemo_v1_update_071:RB: add iom for restart and reorganization of restart
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.6 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 ) CALL lim_dmp( kt )                          ! Ice damping !
220         !                                                           !-------------!
221
222         !                                                           !--------------------!
223         CALL lim_thd( kt )                                          ! Ice thermodynamics !
224         !                                                           !--------------------!
225         IF(ln_ctl) THEN
226            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  4 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
227            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   4 : ', tab2d_2=sist  , clinfo2=' sist      : ')
228         ENDIF
229
230
231         ! Mass and heat fluxes from ice to ocean
232         !                                                           !------------------------------!
233         CALL lim_flx                                                ! Ice/Ocean Mass & Heat fluxes !
234         !                                                           !------------------------------!
235
236         IF( MOD( kt + nfice -1, ninfo ) == 0 .OR. ntmoy == 1 )  THEN        !-----------------!
237            CALL lim_dia( kt )                                    ! Ice Diagnostics !
238         ENDIF                                                       !-----------------!
239
240         !                                                           !-------------!
241         CALL lim_wri( kt )                                          ! Ice outputs !
242         !                                                           !-------------!
243
244         !                                                           !------------------------!
245         IF( lrst_ice ) CALL lim_rst_write( kt )                  ! Write Ice restart file !
246         !                                                           !------------------------!
247
248         ! Re-initialization of forcings
249         qsr_oce (:,:) = 0.e0
250         qsr_ice (:,:) = 0.e0
251         qnsr_oce(:,:) = 0.e0
252         qnsr_ice(:,:) = 0.e0 
253         dqns_ice(:,:) = 0.e0 
254         tprecip (:,:) = 0.e0 
255         sprecip (:,:) = 0.e0
256         fr1_i0  (:,:) = 0.e0
257         fr2_i0  (:,:) = 0.e0
258         evap    (:,:) = 0.e0
259#if defined key_coupled 
260         rrunoff (:,:) = 0.e0
261         calving (:,:) = 0.e0
262#else
263         qla_ice (:,:) = 0.e0
264         dqla_ice(:,:) = 0.e0
265#endif
266      ENDIF
267      !
268   END SUBROUTINE ice_stp
269
270#else
271   !!----------------------------------------------------------------------
272   !!   Default option           Dummy module          NO LIM sea-ice model
273   !!----------------------------------------------------------------------
274CONTAINS
275   SUBROUTINE ice_stp ( kt )     ! Dummy routine
276      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
277   END SUBROUTINE ice_stp
278#endif
279
280   !!======================================================================
281END MODULE icestp
Note: See TracBrowser for help on using the repository browser.