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