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

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

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