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

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

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