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.
icestp1d.F90 in branches/dev_001_SBC/NEMO/C1D_SRC – NEMO

source: branches/dev_001_SBC/NEMO/C1D_SRC/icestp1d.F90 @ 881

Last change on this file since 881 was 881, checked in by ctlod, 16 years ago

dev_001_SBC: Step I: change cpp ket name key_ice_lim into key_lim2 & change names inside modules with extension _2, see ticket: #110

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 16.1 KB
Line 
1MODULE icestp1d
2   !!======================================================================
3   !!                       ***  MODULE icestp1d   ***
4   !!   Sea-Ice model : 1D LIM Sea ice model time-stepping
5   !!======================================================================
6   !! History :   9.0  !  04-10  (C. Ethe)  from icestp, 1D configuration
7   !!----------------------------------------------------------------------
8#if defined key_cfg_1d && defined key_lim2
9   !!----------------------------------------------------------------------
10   !!   'key_cfg_1d'  .AND.                                1D Configuration
11   !!   'key_lim2'                                        Lim sea-ice model
12   !!----------------------------------------------------------------------
13   !!----------------------------------------------------------------------
14   !!   ice_stp_1d       : sea-ice model time-stepping
15   !!----------------------------------------------------------------------
16   USE dom_oce         ! ocean space and time domain
17   USE oce             ! dynamics and tracers variables
18   USE in_out_manager  ! I/O manager
19   USE ice_oce         ! ice variables
20   USE dom_ice_2       ! LIM sea-ice domain
21   USE sbc_oce         ! surface boundary condition: ocean
22   USE sbc_ice         ! surface boundary condition: ice
23   USE daymod          ! calendar
24   USE phycst          ! Define parameters for the routines
25   USE ice_2           ! ice variables
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27   USE limthd_2
28   USE limwri_2
29   USE limrst_2
30
31   USE tradmp          ! damping salinity trend
32   USE dtatem          ! ocean temperature data
33   USE dtasal          ! ocean salinity data
34   USE ocfzpt          ! surface ocean freezing point
35   USE prtctl          ! Print control
36
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   ice_stp_1d  ! called by step.F90
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !!   LIM 2.0 , UCL-LOCEAN-IPSL (2006)
48   !! $Id$
49   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51
52CONTAINS
53
54   SUBROUTINE ice_stp_1d ( kt )
55      !!---------------------------------------------------------------------
56      !!                  ***  ROUTINE ice_stp_1d  ***
57      !!                   
58      !! ** Purpose :   Louvain la Neuve Sea Ice Model time stepping
59      !!
60      !! ** Action  : - call the ice dynamics routine
61      !!              - call the ice advection/diffusion routine
62      !!              - call the ice thermodynamics routine
63      !!              - call the routine that computes mass and
64      !!                heat fluxes at the ice/ocean interface
65      !!              - save the outputs
66      !!              - save the outputs for restart when necessary
67      !!----------------------------------------------------------------------
68      INTEGER, INTENT(in) ::   kt         ! ocean time-step index
69
70      INTEGER                      ::   ji, jj                        ! dummy loop indices
71      REAL(wp)                     ::   ztair2                        ! temporary scalar
72      REAL(wp), DIMENSION(jpi,jpj) ::   zsss_io, zsss2_io, zsss3_io   ! tempory workspaces
73      !!----------------------------------------------------------------------
74
75      IF( kt == nit000 ) THEN
76         IF(lwp) WRITE(numout,*)
77         IF(lwp) WRITE(numout,*) 'ice_stp_1d : Louvain la Neuve Ice Model (LIM)' 
78         IF(lwp) WRITE(numout,*) '~~~~~~~   forced case using bulk formulea'
79         !  Initialize fluxes fields
80         gtaux(:,:) = 0.e0
81         gtauy(:,:) = 0.e0
82      ENDIF
83
84      ! Temperature , salinity and horizonta wind
85      ! sst_io  and sss_io, u_io and v_io  are initialized at nit000 in limistate.F90 (or limrst.F90) with :
86      !              sst_io = sst_io + (nfice - 1) * (tn(:,:,1)+rt0 )
87      !              sss_io = sss_io + (nfice - 1) * sn(:,:,1)
88      !              u_io   = u_io   + (nfice - 1) * un(:,:,1)
89      !              v_io   = v_io   + (nfice - 1) * vn(:,:,1)
90      !    cumulate fields
91      !
92      sst_io(:,:) = sst_io(:,:) + tn(:,:,1) + rt0
93      sss_io(:,:) = sss_io(:,:) + sn(:,:,1)
94      u_io  (:,:) = u_io  (:,:) + un(:,:,1)
95      v_io  (:,:) = v_io  (:,:) + vn(:,:,1)
96 
97     
98      IF( MOD( kt-1, nfice ) == 0 ) THEN
99         
100         ! The LIM model is going to be call
101         sst_io(:,:) = sst_io(:,:) / FLOAT( nfice ) * tmask(:,:,1)
102         sss_io(:,:) = sss_io(:,:) / FLOAT( nfice )
103         u_io  (:,:) = u_io  (:,:) / FLOAT( nfice )
104         v_io  (:,:) = v_io  (:,:) / FLOAT( nfice )
105         gtaux (:,:) = utau  (:,:)
106         gtauy (:,:) = vtau  (:,:)
107
108         zsss_io (:,:) = SQRT( sss_io(:,:) ) 
109         zsss2_io(:,:) =  sss_io(:,:) *  sss_io(:,:)
110         zsss3_io(:,:) = zsss_io(:,:) * zsss_io(:,:) * zsss_io(:,:)
111
112         DO jj = 1, jpj
113            DO ji = 1, jpi
114               tfu(ji,jj)  = ABS ( rt0 - 0.0575       *   sss_io(ji,jj)   &
115                  &                    + 1.710523e-03 * zsss3_io(ji,jj)   &
116                  &                    - 2.154996e-04 * zsss2_io(ji,jj) ) * tms(ji,jj)
117            END DO
118         END DO
119         
120         
121         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
122            CALL prt_ctl_info('Ice Forcings ')
123            CALL prt_ctl(tab2d_1=qsr_oce ,clinfo1=' qsr_oce  : ', tab2d_2=qsr_ice , clinfo2=' qsr_ice   : ')
124            CALL prt_ctl(tab2d_1=qnsr_oce,clinfo1=' qnsr_oce : ', tab2d_2=qnsr_ice, clinfo2=' qnsr_ice  : ')
125            CALL prt_ctl(tab2d_1=evap    ,clinfo1=' evap     : ')
126            CALL prt_ctl(tab2d_1=tprecip ,clinfo1=' precip   : ', tab2d_2=sprecip , clinfo2=' Snow      : ')
127            CALL prt_ctl(tab2d_1=gtaux   ,clinfo1=' u-stress : ', tab2d_2=gtauy   , clinfo2=' v-stress  : ')
128            CALL prt_ctl(tab2d_1=sst_io  ,clinfo1=' sst      : ', tab2d_2=sss_io  , clinfo2=' sss       : ')
129            CALL prt_ctl(tab2d_1=u_io    ,clinfo1=' u_io     : ', tab2d_2=v_io    , clinfo2=' v_io      : ')
130            CALL prt_ctl(tab2d_1=hsnif   ,clinfo1=' hsnif  1 : ', tab2d_2=hicif   , clinfo2=' hicif     : ')
131            CALL prt_ctl(tab2d_1=frld    ,clinfo1=' frld   1 : ', tab2d_2=sist    , clinfo2=' sist      : ')
132         ENDIF
133                                                                     
134         DO jj = 1, jpj
135            DO ji = 1, jpi
136               tio_u(ji,jj) = - gtaux(ji,jj) / rau0
137               tio_v(ji,jj) = - gtauy(ji,jj) / rau0               
138               ztair2       = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj)           
139               ust2s(ji,jj) = ( SQRT( ztair2  )  / rau0 )  * tms(ji,jj) 
140            END DO
141         END DO
142
143         !                                                           !-----------------------!
144         CALL lim_rst_opn( kt )                                      ! Open Ice restart file !
145         !                                                           !-----------------------!
146
147         !                                                           !--------------------!
148         CALL lim_thd( kt )                                          ! Ice thermodynamics !
149         !                                                           !--------------------!
150         IF(ln_ctl) THEN
151            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  2 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
152            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   2 : ', tab2d_2=sist  , clinfo2=' sist      : ')
153            CALL prt_ctl(tab2d_1=u_io  ,clinfo1=' u_io   4 : ', tab2d_2=v_io  , clinfo2=' v_io      : ')
154            CALL prt_ctl(tab2d_1=tio_u  ,clinfo1=' tio_u  4 : ', tab2d_2=tio_v  , clinfo2=' tio_v     : ')
155         ENDIF
156
157
158
159         ! Mass and heat fluxes from ice to ocean
160         !                                                           !------------------------------!
161         CALL lim_flx                                                ! Ice/Ocean Mass & Heat fluxes !
162         !                                                           !------------------------------!
163
164         IF(ln_ctl) THEN
165            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  7 : ', tab2d_2=hicif , clinfo2=' hicif   : ')
166            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   7 : ', tab2d_2=sist  , clinfo2=' sist      : ')
167            CALL prt_ctl(tab2d_1=tio_u  ,clinfo1=' tio_u  7 : ', tab2d_2=tio_v  , clinfo2=' tio_v     : ')
168         ENDIF
169         !                                                           !-------------!
170         CALL lim_wri( kt )                                          ! Ice outputs !
171         !                                                           !-------------!
172
173         !                                                           !------------------------!
174         IF( lrst_ice ) CALL lim_rst_write( kt )                     ! Write Ice restart file !
175         !                                                           !------------------------!
176
177         ! Re-initialization of forcings
178         qsr_oce (:,:) = 0.e0
179         qsr_ice (:,:) = 0.e0
180         qnsr_oce(:,:) = 0.e0
181         qnsr_ice(:,:) = 0.e0 
182         dqns_ice(:,:) = 0.e0 
183         tprecip (:,:) = 0.e0 
184         sprecip (:,:) = 0.e0
185         qla_ice (:,:) = 0.e0
186         dqla_ice(:,:) = 0.e0
187         fr1_i0  (:,:) = 0.e0
188         fr2_i0  (:,:) = 0.e0
189         evap    (:,:) = 0.e0
190
191        CALL oce_sbc_1d ( kt )
192
193      ENDIF
194
195   END SUBROUTINE ice_stp_1d
196   
197   
198   SUBROUTINE oce_sbc_1d( kt )
199      !!---------------------------------------------------------------------
200      !!                   ***  ROUTINE oce_sbc_1d  ***
201      !!                   
202      !! ** Purpose : - Ocean surface boundary conditions with LIM sea-ice
203      !!        model in forced mode using bulk formulea
204      !!----------------------------------------------------------------------
205      INTEGER, INTENT(in) ::   kt   ! ocean time step
206      !
207      INTEGER  ::   ji, jj                   ! dummy loop indices
208      REAL(wp) ::   ztxy
209      !!----------------------------------------------------------------------
210
211      ! 1. initialization to zero at kt = nit000
212      ! ---------------------------------------
213     
214      IF( kt == nit000 ) THEN     
215         qsr    (:,:) = 0.e0
216         qns    (:,:) = 0.e0
217         qrp    (:,:) = 0.e0
218         emp    (:,:) = 0.e0
219         emps   (:,:) = 0.e0
220         erp    (:,:) = 0.e0
221#if ! defined key_dynspg_rl
222         dmp    (:,:) = 0.e0
223#endif
224      ENDIF
225
226      CALL oce_sbc_dmp       ! Computation of internal and evaporation damping terms       
227
228      ! Surface Ocean fluxes
229      ! ====================
230     
231      ! Surface heat flux (W/m2)
232      ! -----------------
233     
234      qns (:,:) = fnsolar(:,:)                    ! non solar heat flux
235      qsr (:,:) = fsolar(:,:)                     ! solar flux
236     
237#if ! defined key_dynspg_rl     
238      ! total concentration/dilution effect (use on SSS)
239      emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:) + empold
240     
241      ! total volume flux (use on sea-surface height)
242      emp (:,:) = fmass(:,:) -   dmp(:,:) + runoff(:,:) + erp(:,:) + empold     
243#else
244      ! Rigid-lid (emp=emps=E-P-R+Erp)
245      emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)     ! freshwater flux
246      emp (:,:) = emps(:,:)
247     
248#endif
249     
250      ! Surface stress
251      ! --------------
252     
253      ! update the stress beloww sea-ice area
254      DO jj = 1, jpjm1
255         DO ji = 1, fs_jpim1   ! vertor opt.
256            ztxy        = freezn(ji,jj)             ! ice/ocean indicator at T-points
257            utau(ji,jj) = (1.-ztxy) * utau(ji,jj) + ztxy * ftaux(ji,jj)    ! stress at the ocean surface
258            vtau(ji,jj) = (1.-ztxy) * vtau(ji,jj) + ztxy * ftauy(ji,jj)
259         END DO
260      END DO
261     
262      ! boundary condition on the stress (utau,vtau)
263      CALL lbc_lnk( utau, 'U', -1. )
264      CALL lbc_lnk( vtau, 'V', -1. )
265     
266      ! Re-initialization of fluxes
267      sst_io(:,:) = 0.e0
268      sss_io(:,:) = 0.e0
269      u_io  (:,:) = 0.e0
270      v_io  (:,:) = 0.e0
271      !
272   END SUBROUTINE oce_sbc_1d
273   
274#if defined key_dtasal
275   !!----------------------------------------------------------------------
276   !!   'key_dtasal'                                          salinity data
277   !!----------------------------------------------------------------------
278   SUBROUTINE oce_sbc_dmp
279      !!---------------------------------------------------------------------
280      !!                   ***  ROUTINE oce_sbc_dmp  ***
281      !!                   
282      !! ** Purpose : Computation of internal and evaporation damping terms
283      !!        for ocean surface boundary conditions
284      !!----------------------------------------------------------------------
285      INTEGER  ::   ji, jj                   ! dummy loop indices
286      REAL(wp) ::   zerp, zsrp
287#if ! defined key_dynspg_rl
288      REAL(wp) ::   zwei
289      REAL(wp) ::   zerpplus(jpi,jpj), zerpminus(jpi,jpj)
290      REAL(wp) ::   zplus, zminus, zadefi
291# if defined key_tradmp
292      INTEGER jk
293      REAL(wp), DIMENSION(jpi,jpj) ::   zstrdmp
294# endif
295#endif
296      !!----------------------------------------------------------------------
297
298      ! sea ice indicator (1 or 0)
299      DO jj = 1, jpj
300         DO ji = 1, jpi
301            freezn(ji,jj) = MAX(0., SIGN(1., freeze(ji,jj)-rsmall) )
302         END DO
303      END DO
304
305      ! Initialisation
306      ! --------------
307      ! Restoring coefficients on SST and SSS   
308      zsrp = dqdt0 * ro0cpr * rauw   ! (Kg/m2/s)
309
310#if ! defined key_dynspg_rl
311      ! Free-surface
312         
313      ! Internal damping
314# if defined key_tradmp
315      ! Vertical mean of dampind trend (computed in tradmp module)
316      zstrdmp(:,:) = 0.e0
317      DO jk = 1, jpk
318         zstrdmp(:,:) = zstrdmp(:,:) + strdmp(:,:,jk) * fse3t(:,:,jk)
319      END DO
320      ! volume flux associated to internal damping to climatology
321      dmp(:,:) = zstrdmp(:,:) * rauw / ( sss_io(:,:) + 1.e-20 )
322# else
323      dmp(:,:) = 0.e0            ! No internal damping
324# endif
325     
326      !   evaporation damping term ( Surface restoring )
327      zerpplus (:,:) = 0.e0
328      zerpminus(:,:) = 0.e0
329      zplus          =  15. / rday
330      zminus         = -15. / rday
331     
332      DO jj = 1, jpj
333         DO ji = 1, jpi
334            zerp = ( 1. - 2.*upsrnfh(ji,jj) ) * zsrp   &
335               & * ( sss_io(ji,jj) - s_dta(ji,jj,1) )     &
336               & / ( sss_io(ji,jj) + 1.e-20        )
337            erp(ji,jj) = zerp
338            zerpplus (ji,jj) = MAX( erp(ji,jj), 0.e0 )
339            zerpminus(ji,jj) = MIN( erp(ji,jj), 0.e0 )
340         END DO
341      END DO
342
343      aplus  = 0.e0
344      aminus = 0.e0
345      DO jj = 1, jpj
346         DO ji = 1, jpi
347            zwei   = e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
348            aplus  = aplus  + zerpplus (ji,jj) * zwei
349            aminus = aminus - zerpminus(ji,jj) * zwei
350         END DO
351      END DO
352
353      IF(ln_ctl)   WRITE(numout,*) ' oce_sbc_dmp : a+ = ', aplus, ' a- = ', aminus
354#else
355      ! Rigid-lid (emp=emps=E-P-R+Erp)
356     
357      erp(:,:) = ( 1. - freezn(:,:) ) * zsrp    &   ! surface restoring term
358         &     * ( sss_io(:,:) - s_dta(:,:,1) )     &
359         &     / ( sss_io(:,:) + 1.e-20      )
360#endif
361
362   END SUBROUTINE oce_sbc_dmp
363   
364#else
365   !!----------------------------------------------------------------------
366   !!   Dummy routine                                      NO salinity data
367   !!----------------------------------------------------------------------
368   SUBROUTINE oce_sbc_dmp         ! Dummy routine
369      WRITE(*,*) 'oce_sbc_dmp: you should not have seen that print! error?'
370   END SUBROUTINE oce_sbc_dmp
371#endif
372#else
373   !!----------------------------------------------------------------------
374   !!   Default option           Dummy module   NO 1D && NO LIM sea-ice model
375   !!----------------------------------------------------------------------
376CONTAINS
377   SUBROUTINE ice_stp_1d ( kt )     ! Dummy routine
378      WRITE(*,*) 'ice_stp_1d: You should not have seen this print! error?', kt
379   END SUBROUTINE ice_stp_1d
380   SUBROUTINE oce_sbc_1d ( kt )     ! Dummy routine
381      WRITE(*,*) 'oce_sbc_1d: You should not have seen this print! error?', kt
382   END SUBROUTINE oce_sbc_1d
383#endif   
384
385   !!======================================================================
386END MODULE icestp1d
Note: See TracBrowser for help on using the repository browser.