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

source: tags/nemo_v1_08/NEMO/C1D_SRC/icestp1d.F90 @ 2143

Last change on this file since 2143 was 321, checked in by opalod, 19 years ago

nemo_v1_bugfix_014:RB: update control print for 1D model

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.3 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   USE prtctl          ! Print control
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(ln_ctl) THEN         ! print mean trends (used for debugging)
141            CALL prt_ctl_info('Ice Forcings ')
142            CALL prt_ctl(tab2d_1=qsr_oce ,clinfo1=' qsr_oce  : ', tab2d_2=qsr_ice , clinfo2=' qsr_ice   : ')
143            CALL prt_ctl(tab2d_1=qnsr_oce,clinfo1=' qnsr_oce : ', tab2d_2=qnsr_ice, clinfo2=' qnsr_ice  : ')
144            CALL prt_ctl(tab2d_1=evap    ,clinfo1=' evap     : ')
145            CALL prt_ctl(tab2d_1=tprecip ,clinfo1=' precip   : ', tab2d_2=sprecip , clinfo2=' Snow      : ')
146            CALL prt_ctl(tab2d_1=gtaux   ,clinfo1=' u-stress : ', tab2d_2=gtauy   , clinfo2=' v-stress  : ')
147            CALL prt_ctl(tab2d_1=sst_io  ,clinfo1=' sst      : ', tab2d_2=sss_io  , clinfo2=' sss       : ')
148            CALL prt_ctl(tab2d_1=u_io    ,clinfo1=' u_io     : ', tab2d_2=v_io    , clinfo2=' v_io      : ')
149            CALL prt_ctl(tab2d_1=hsnif   ,clinfo1=' hsnif  1 : ', tab2d_2=hicif   , clinfo2=' hicif     : ')
150            CALL prt_ctl(tab2d_1=frld    ,clinfo1=' frld   1 : ', tab2d_2=sist    , clinfo2=' sist      : ')
151         ENDIF
152
153         ! Ice model call
154         numit = numit + nfice                                       ! Friction velocity
155                                                                     
156         DO jj = 1, jpj
157            DO ji = 1, jpi
158               tio_u(ji,jj) = - gtaux(ji,jj) / rau0
159               tio_v(ji,jj) = - gtauy(ji,jj) / rau0               
160               ztair2       = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj)           
161               ust2s(ji,jj) = ( SQRT( ztair2  )  / rau0 )  * tms(ji,jj) 
162            END DO
163         END DO
164
165         !                                                           !--------------------!
166         CALL lim_thd                                                ! Ice thermodynamics !
167         !                                                           !--------------------!
168         IF(ln_ctl) THEN
169            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  2 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
170            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   2 : ', tab2d_2=sist  , clinfo2=' sist      : ')
171            CALL prt_ctl(tab2d_1=u_io  ,clinfo1=' u_io   4 : ', tab2d_2=v_io  , clinfo2=' v_io      : ')
172            CALL prt_ctl(tab2d_1=tio_u  ,clinfo1=' tio_u  4 : ', tab2d_2=tio_v  , clinfo2=' tio_v     : ')
173         ENDIF
174
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(ln_ctl) THEN
183            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  7 : ', tab2d_2=hicif , clinfo2=' hicif   : ')
184            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   7 : ', tab2d_2=sist  , clinfo2=' sist      : ')
185            CALL prt_ctl(tab2d_1=tio_u  ,clinfo1=' tio_u  7 : ', tab2d_2=tio_v  , clinfo2=' tio_v     : ')
186         ENDIF
187         !                                                           !-------------!
188         CALL lim_wri                                                ! Ice outputs !
189         !                                                           !-------------!
190
191         IF( MOD( numit, nstock ) == 0 .OR. numit == nlast ) THEN
192            !                                                        !------------------!
193            CALL lim_rst_write( numit )                              ! Ice restart file !
194            !                                                        !------------------!
195         ENDIF
196
197         ! Re-initialization of forcings
198         qsr_oce (:,:) = 0.e0
199         qsr_ice (:,:) = 0.e0
200         qnsr_oce(:,:) = 0.e0
201         qnsr_ice(:,:) = 0.e0 
202         dqns_ice(:,:) = 0.e0 
203         tprecip (:,:) = 0.e0 
204         sprecip (:,:) = 0.e0
205         qla_ice (:,:) = 0.e0
206         dqla_ice(:,:) = 0.e0
207         fr1_i0  (:,:) = 0.e0
208         fr2_i0  (:,:) = 0.e0
209         evap    (:,:) = 0.e0
210
211        CALL oce_sbc_1d ( kt )
212
213      ENDIF
214
215   END SUBROUTINE ice_stp_1d
216   
217   
218   SUBROUTINE oce_sbc_1d( kt )
219      !!---------------------------------------------------------------------
220      !!                   ***  ROUTINE oce_sbc_1d  ***
221      !!                   
222      !! ** Purpose : - Ocean surface boundary conditions with LIM sea-ice
223      !!        model in forced mode using bulk formulea
224      !!
225      !! History :
226      !!   1.0  !  99-11  (M. Imbard)  Original code
227      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
228      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
229      !!----------------------------------------------------------------------
230      !! * arguments
231      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
232
233      !! * Local declarations
234      INTEGER  ::   ji, jj                   ! dummy loop indices
235      REAL(wp) ::   ztxy
236      !!----------------------------------------------------------------------
237
238      ! 1. initialization to zero at kt = nit000
239      ! ---------------------------------------
240     
241      IF( kt == nit000 ) THEN     
242         qsr    (:,:) = 0.e0
243         qt     (:,:) = 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      qt  (:,:) = fnsolar(:,:) + fsolar(:,:)     ! non solar heat flux + solar flux
262      qsr (:,:) = fsolar(:,:)                     ! solar flux
263     
264#if defined key_dynspg_fsc     
265      ! total concentration/dilution effect (use on SSS)
266      emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:) + empold
267     
268      ! total volume flux (use on sea-surface height)
269      emp (:,:) = fmass(:,:) -   dmp(:,:) + runoff(:,:) + erp(:,:) + empold     
270#else
271      ! Rigid-lid (emp=emps=E-P-R+Erp)
272      emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)     ! freshwater flux
273      emp (:,:) = emps(:,:)
274     
275#endif
276     
277      ! Surface stress
278      ! --------------
279     
280      ! update the stress beloww sea-ice area
281      DO jj = 1, jpjm1
282         DO ji = 1, fs_jpim1   ! vertor opt.
283            ztxy        = freezn(ji,jj)             ! ice/ocean indicator at T-points
284            taux(ji,jj) = (1.-ztxy) * taux(ji,jj) + ztxy * ftaux(ji,jj)    ! stress at the ocean surface
285            tauy(ji,jj) = (1.-ztxy) * tauy(ji,jj) + ztxy * ftauy(ji,jj)
286         END DO
287      END DO
288     
289      ! boundary condition on the stress (taux,tauy)
290      CALL lbc_lnk( taux, 'U', -1. )
291      CALL lbc_lnk( tauy, 'V', -1. )
292     
293      ! Re-initialization of fluxes
294      sst_io(:,:) = 0.e0
295      sss_io(:,:) = 0.e0
296      u_io  (:,:) = 0.e0
297      v_io  (:,:) = 0.e0
298     
299     
300   END SUBROUTINE oce_sbc_1d
301   
302#if defined key_dtasal
303   !!----------------------------------------------------------------------
304   !!   'key_dtasal'                                          salinity data
305   !!----------------------------------------------------------------------
306   SUBROUTINE oce_sbc_dmp
307      !!---------------------------------------------------------------------
308      !!                   ***  ROUTINE oce_sbc_dmp  ***
309      !!                   
310      !! ** Purpose : Computation of internal and evaporation damping terms
311      !!        for ocean surface boundary conditions
312      !!
313      !! History :
314      !!   9.0  !  04-01  (G. Madec, C. Ethe)  Original code
315      !!----------------------------------------------------------------------
316      !! * Local declarations
317      INTEGER ::   ji, jj                   ! dummy loop indices
318      REAL(wp) ::   zerp, ztrp, zsrp
319#if defined key_dynspg_fsc
320      REAL(wp) ::   zwei
321      REAL(wp) ::   zerpplus(jpi,jpj), zerpminus(jpi,jpj)
322      REAL(wp) ::   zplus, zminus, zadefi
323# if defined key_tradmp
324      INTEGER jk
325      REAL(wp), DIMENSION(jpi,jpj) ::   zstrdmp
326# endif
327#endif
328      !!----------------------------------------------------------------------
329
330
331      ! sea ice indicator (1 or 0)
332      DO jj = 1, jpj
333         DO ji = 1, jpi
334            freezn(ji,jj) = MAX(0., SIGN(1., freeze(ji,jj)-rsmall) )
335         END DO
336      END DO
337
338      ! Initialisation
339      ! --------------
340      ! Restoring coefficients on SST and SSS   
341      ztrp  = -40.                   ! (W/m2/K)   
342      zsrp  = ztrp * ro0cpr * rauw   ! (Kg/m2/s2)
343
344#if defined key_dynspg_fsc 
345      ! Free-surface
346         
347      ! Internal damping
348# if defined key_tradmp
349      ! Vertical mean of dampind trend (computed in tradmp module)
350      zstrdmp(:,:) = 0.e0
351      DO jk = 1, jpk
352         zstrdmp(:,:) = zstrdmp(:,:) + strdmp(:,:,jk) * fse3t(:,:,jk)
353      END DO
354      ! volume flux associated to internal damping to climatology
355      dmp(:,:) = zstrdmp(:,:) * rauw / ( sss_io(:,:) + 1.e-20 )
356# else
357      dmp(:,:) = 0.e0            ! No internal damping
358# endif
359     
360      !   evaporation damping term ( Surface restoring )
361      zerpplus (:,:) = 0.e0
362      zerpminus(:,:) = 0.e0
363      zplus          =  15. / rday
364      zminus         = -15. / rday
365     
366      DO jj = 1, jpj
367         DO ji = 1, jpi
368            zerp = ( 1. - 2.*upsrnfh(ji,jj) ) * zsrp   &
369               & * ( sss_io(ji,jj) - s_dta(ji,jj,1) )     &
370               & / ( sss_io(ji,jj) + 1.e-20        )
371            erp(ji,jj) = zerp
372            zerpplus (ji,jj) = MAX( erp(ji,jj), 0.e0 )
373            zerpminus(ji,jj) = MIN( erp(ji,jj), 0.e0 )
374         END DO
375      END DO
376
377      aplus  = 0.e0
378      aminus = 0.e0
379      DO jj = 1, jpj
380         DO ji = 1, jpi
381            zwei   = e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
382            aplus  = aplus  + zerpplus (ji,jj) * zwei
383            aminus = aminus - zerpminus(ji,jj) * zwei
384         END DO
385      END DO
386
387      IF(ln_ctl)   WRITE(numout,*) ' oce_sbc_dmp : a+ = ', aplus, ' a- = ', aminus
388#else
389      ! Rigid-lid (emp=emps=E-P-R+Erp)
390     
391      erp(:,:) = ( 1. - freezn(:,:) ) * zsrp    &   ! surface restoring term
392         &     * ( sss_io(:,:) - s_dta(:,:,1) )     &
393         &     / ( sss_io(:,:) + 1.e-20      )
394#endif
395
396   END SUBROUTINE oce_sbc_dmp
397   
398#else
399   !!----------------------------------------------------------------------
400   !!   Dummy routine                                      NO salinity data
401   !!----------------------------------------------------------------------
402   SUBROUTINE oce_sbc_dmp         ! Dummy routine
403      WRITE(*,*) 'oce_sbc_dmp: you should not have seen that print! error?'
404   END SUBROUTINE oce_sbc_dmp
405#endif
406#else
407   !!----------------------------------------------------------------------
408   !!   Default option           Dummy module   NO 1D && NO LIM sea-ice model
409   !!----------------------------------------------------------------------
410CONTAINS
411   SUBROUTINE ice_stp_1d ( kt )     ! Dummy routine
412      WRITE(*,*) 'ice_stp_1d: You should not have seen this print! error?', kt
413   END SUBROUTINE ice_stp_1d
414   SUBROUTINE oce_sbc_1d ( kt )     ! Dummy routine
415      WRITE(*,*) 'oce_sbc_1d: You should not have seen this print! error?', kt
416   END SUBROUTINE oce_sbc_1d
417#endif   
418   !!======================================================================
419END MODULE icestp1d
Note: See TracBrowser for help on using the repository browser.