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

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

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