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

source: tags/nemo_v1_02/NEMO/C1D_SRC/icestp1d.F90 @ 7041

Last change on this file since 7041 was 253, checked in by opalod, 19 years ago

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