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 on Ticket #255 – Attachment – NEMO

Ticket #255: icestp.F90

File icestp.F90, 13.5 KB (added by nemo_user, 16 years ago)
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#if defined key_agrif
40   USE agrif_opa_interp
41   USE agrif_opa_update
42#endif
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   ice_stp    ! called by step.F90
48
49   !! * Substitutions
50#  include "domzgr_substitute.h90"
51#  include "vectopt_loop_substitute.h90"
52   !!-----------------------------------------------------
53   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)
54   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/icestp.F90,v 1.9 2006/11/15 10:45:39 opalod Exp $
55   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
56   !!-----------------------------------------------------
57
58CONTAINS
59
60   SUBROUTINE ice_stp ( kt )
61      !!---------------------------------------------------------------------
62      !!                  ***  ROUTINE ice_stp  ***
63      !!                   
64      !! ** Purpose :   Louvain la Neuve Sea Ice Model time stepping
65      !!
66      !! ** Action  : - call the ice dynamics routine
67      !!              - call the ice advection/diffusion routine
68      !!              - call the ice thermodynamics routine
69      !!              - call the routine that computes mass and
70      !!                heat fluxes at the ice/ocean interface
71      !!              - save the outputs
72      !!              - save the outputs for restart when necessary
73      !!----------------------------------------------------------------------
74      INTEGER, INTENT(in)          ::   kt       ! ocean time-step index
75
76      INTEGER                      ::   ji, jj   ! dummy loop indices
77      REAL(wp), DIMENSION(jpi,jpj) ::   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! FD debug ice-ocean stress
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            u_io(ji,jj) = u_io(ji,jj) + 0.5 * ( ub(ji-1,jj  ,1) + ub(ji-1,jj-1,1) )
114            v_io(ji,jj) = v_io(ji,jj) + 0.5 * ( vb(ji  ,jj-1,1) + vb(ji-1,jj-1,1) )
115         END DO
116      END DO
117     
118     
119      IF( MOD( kt-1, nfice ) == 0 ) THEN
120         
121         ! The LIM model is going to be call
122         sst_io(:,:) = sst_io(:,:) / FLOAT( nfice ) * tmask(:,:,1)
123         sss_io(:,:) = sss_io(:,:) / FLOAT( nfice )
124
125         ! stress from ocean U- and V-points to ice U,V point
126         DO jj = 2, jpj
127            DO ji = fs_2, jpi   ! vector opt.
128!Sujie               gtaux(ji,jj) = 0.5 * ( taux(ji-1,jj  ) + taux(ji-1,jj-1) )
129!Sujie               gtauy(ji,jj) = 0.5 * ( tauy(ji  ,jj-1) + tauy(ji-1,jj-1) )
130! FD modifications from JM Molines
131               gtaux(ji,jj) = 0.5 * ( tauxwi(ji-1,jj  ) + tauxwi(ji-1,jj-1) )
132               gtauy(ji,jj) = 0.5 * ( tauywi(ji  ,jj-1) + tauywi(ji-1,jj-1) )
133               u_io  (ji,jj) = u_io(ji,jj) / FLOAT( nfice )
134               v_io  (ji,jj) = v_io(ji,jj) / FLOAT( nfice )
135            END DO
136         END DO
137
138         ! lateral boundary condition
139         CALL lbc_lnk( gtaux(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
140         CALL lbc_lnk( gtauy(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
141         CALL lbc_lnk( u_io (:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
142         CALL lbc_lnk( v_io (:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
143
144!!gmbug  in the ocean freezing point computed as :
145!!gm           fzptn (ji,jj) = ( -0.0575 + 1.710523e-3 * SQRT( sn(ji,jj,1) )   &
146!!gm                                     - 2.154996e-4 *       sn(ji,jj,1)   ) * sn(ji,jj,1)   !!   &
147!!gm           !!                        - 7.53e-4 * pressure
148!!gm
149!!gm!bug this is much more accurate and efficient computation
150!!gm       **************************************************
151!!gm freezing point from unesco:
152!!gm     real function tf(s,p)
153!!gm   function to compute the freezing point of seawater
154!!gm
155!!gm   reference: unesco tech. papers in the marine science no. 28. 1978
156!!gm   eighth report jpots
157!!gm   annex 6 freezing point of seawater f.j. millero pp.29-35.
158!!gm
159!!gm  units:
160!!gm         pressure      p          decibars
161!!gm         salinity      s          pss-78
162!!gm         temperature   tf         degrees celsius
163!!gm         freezing pt.
164!!gm************************************************************
165!!gm  checkvalue: tf= -2.588567 deg. c for s=40.0, p=500. decibars
166!!gm     tf=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p
167!!gm     return
168!!gm     end
169!!gm!bug
170
171
172!!gm      DO jj = 1, jpj
173!!gm         DO ji = 1, jpi
174!!gm            tfu(ji,jj)  = (  rt0 + ( - 0.0575                              &
175!!gm               &                     + 1.710523e-3 * SQRT( sss_io(ji,jj) )   &
176!!gm               &                     - 2.154996e-4 *       sss_io(ji,jj)   ) * sss_io(ji,jj)  ) * tms(ji,jj)
177!!gm         END DO
178!!gm      END DO
179!!gm
180      zsss_io (:,:) = SQRT( sss_io(:,:) ) 
181      zsss2_io(:,:) =  sss_io(:,:) *  sss_io(:,:)
182      zsss3_io(:,:) = zsss_io(:,:) * zsss_io(:,:) * zsss_io(:,:)
183
184      DO jj = 1, jpj
185         DO ji = 1, jpi
186            tfu(ji,jj)  = ABS ( rt0 - 0.0575       *   sss_io(ji,jj)   &
187               &                    + 1.710523e-03 * zsss3_io(ji,jj)   &
188               &                    - 2.154996e-04 * zsss2_io(ji,jj) ) * tms(ji,jj)
189         END DO
190      END DO
191   
192         
193
194         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
195            CALL prt_ctl_info('Ice Forcings ')
196            CALL prt_ctl(tab2d_1=qsr_oce ,clinfo1=' qsr_oce  : ', tab2d_2=qsr_ice , clinfo2=' qsr_ice   : ')
197            CALL prt_ctl(tab2d_1=qnsr_oce,clinfo1=' qnsr_oce : ', tab2d_2=qnsr_ice, clinfo2=' qnsr_ice  : ')
198            CALL prt_ctl(tab2d_1=evap    ,clinfo1=' evap     : ')
199            CALL prt_ctl(tab2d_1=tprecip ,clinfo1=' precip   : ', tab2d_2=sprecip , clinfo2=' Snow      : ')
200            CALL prt_ctl(tab2d_1=gtaux   ,clinfo1=' u-stress : ', tab2d_2=gtauy   , clinfo2=' v-stress  : ')
201            CALL prt_ctl(tab2d_1=sst_io  ,clinfo1=' sst      : ', tab2d_2=sss_io  , clinfo2=' sss       : ')
202            CALL prt_ctl(tab2d_1=u_io    ,clinfo1=' u_io     : ', tab2d_2=v_io    , clinfo2=' v_io      : ')
203            CALL prt_ctl(tab2d_1=hsnif   ,clinfo1=' hsnif  1 : ', tab2d_2=hicif   , clinfo2=' hicif     : ')
204            CALL prt_ctl(tab2d_1=frld    ,clinfo1=' frld   1 : ', tab2d_2=sist    , clinfo2=' sist      : ')
205         ENDIF
206
207         !                                                           !-----------------------!
208         CALL lim_rst_opn( kt )                                   ! Open Ice restart file !
209         !                                                           !-----------------------!
210
211         !                                                           !--------------!
212         CALL lim_dyn( kt )                                          ! Ice dynamics !   ( rheology/dynamics )
213         !                                                           !--------------!
214         IF(ln_ctl) THEN
215            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  2 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
216            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   2 : ', tab2d_2=sist  , clinfo2=' sist      : ')
217         ENDIF
218
219         !                                                           !---------------!
220         CALL lim_trp( kt )                                          ! Ice transport !  ( Advection/diffusion )
221         !                                                           !---------------!
222         IF(ln_ctl) THEN
223            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  3 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
224            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   3 : ', tab2d_2=sist  , clinfo2=' sist      : ')
225         ENDIF
226
227         !                                                           !-------------!
228         IF( ln_limdmp ) THEN                                        ! Ice damping !
229         !                                                           !-------------!
230#if defined key_agrif
231            IF( Agrif_Root() )   CALL lim_dmp(kt)
232#else
233            CALL lim_dmp(kt)                                         
234#endif
235         ENDIF
236         !                                                           !--------------------!
237         CALL lim_thd( kt )                                          ! Ice thermodynamics !
238         !                                                           !--------------------!
239         IF(ln_ctl) THEN
240            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  4 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
241            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   4 : ', tab2d_2=sist  , clinfo2=' sist      : ')
242         ENDIF
243
244! FD add-on
245#if defined key_agrif
246         CALL Agrif_tra_ice( kt )
247         IF (.NOT.Agrif_Root())    CALL Agrif_Update_ice( kt )
248#endif
249
250         ! Mass and heat fluxes from ice to ocean
251         !                                                           !------------------------------!
252         CALL lim_flx                                                ! Ice/Ocean Mass & Heat fluxes !
253         !                                                           !------------------------------!
254
255         IF( MOD( kt + nfice -1, ninfo ) == 0 .OR. ntmoy == 1 )  THEN        !-----------------!
256            CALL lim_dia( kt )                                    ! Ice Diagnostics !
257         ENDIF                                                       !-----------------!
258
259         !                                                           !-------------!
260         CALL lim_wri( kt )                                          ! Ice outputs !
261         !                                                           !-------------!
262
263         !                                                           !------------------------!
264         IF( lrst_ice ) CALL lim_rst_write( kt )                  ! Write Ice restart file !
265         !                                                           !------------------------!
266
267         ! Re-initialization of forcings
268         qsr_oce (:,:) = 0.e0
269         qsr_ice (:,:) = 0.e0
270         qnsr_oce(:,:) = 0.e0
271         qnsr_ice(:,:) = 0.e0 
272         dqns_ice(:,:) = 0.e0 
273         tprecip (:,:) = 0.e0 
274         sprecip (:,:) = 0.e0
275         fr1_i0  (:,:) = 0.e0
276         fr2_i0  (:,:) = 0.e0
277         evap    (:,:) = 0.e0
278#if defined key_coupled 
279         rrunoff (:,:) = 0.e0
280         calving (:,:) = 0.e0
281#else
282         qla_ice (:,:) = 0.e0
283         dqla_ice(:,:) = 0.e0
284#endif
285      ENDIF
286      !
287   END SUBROUTINE ice_stp
288
289#else
290   !!----------------------------------------------------------------------
291   !!   Default option           Dummy module          NO LIM sea-ice model
292   !!----------------------------------------------------------------------
293CONTAINS
294   SUBROUTINE ice_stp ( kt )     ! Dummy routine
295      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
296   END SUBROUTINE ice_stp
297#endif
298
299   !!======================================================================
300END MODULE icestp