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 trunk/NEMO/C1D_SRC – NEMO

source: trunk/NEMO/C1D_SRC/icestp1d.F90 @ 833

Last change on this file since 833 was 833, checked in by rblod, 16 years ago

Merge branche dev_002_LIM back to trunk ticket #70 and #71

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