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

source: tags/nemo_v1_03/NEMO/C1D_SRC/icestp1d.F90 @ 4100

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

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