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.
sbcice_lim.F90 in branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90 @ 5728

Last change on this file since 5728 was 4333, checked in by clem, 10 years ago

remove remaining bugs in LIM3, so that it can run in both regional and global config

  • Property svn:keywords set to Id
File size: 44.6 KB
Line 
1MODULE sbcice_lim
2   !!======================================================================
3   !!                       ***  MODULE  sbcice_lim  ***
4   !! Surface module :  update the ocean surface boundary condition over ice
5   !!       &           covered area using LIM sea-ice model
6   !! Sea-Ice model  :  LIM-3 Sea ice model time-stepping
7   !!=====================================================================
8   !! History :  2.0  ! 2006-12  (M. Vancoppenolle) Original code
9   !!            3.0  ! 2008-02  (C. Talandier)  Surface module from icestp.F90
10   !!             -   ! 2008-04  (G. Madec)  sltyle and lim_ctl routine
11   !!            3.3  ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
12   !!            3.4  ! 2011-01  (A Porter)  dynamical allocation
13   !!             -   ! 2012-10  (C. Rousset)  add lim_diahsb
14   !!----------------------------------------------------------------------
15#if defined key_lim3
16   !!----------------------------------------------------------------------
17   !!   'key_lim3' :                                  LIM 3.0 sea-ice model
18   !!----------------------------------------------------------------------
19   !!   sbc_ice_lim  : sea-ice model time-stepping and update ocean sbc over ice-covered area
20   !!   lim_ctl       : alerts in case of ice model crash
21   !!   lim_prt_state : ice control print at a given grid point
22   !!----------------------------------------------------------------------
23   USE oce             ! ocean dynamics and tracers
24   USE dom_oce         ! ocean space and time domain
25   USE par_ice         ! sea-ice parameters
26   USE ice             ! LIM-3: ice variables
27   USE iceini          ! LIM-3: ice initialisation
28   USE dom_ice         ! LIM-3: ice domain
29
30   USE sbc_oce         ! Surface boundary condition: ocean fields
31   USE sbc_ice         ! Surface boundary condition: ice   fields
32   USE sbcblk_core     ! Surface boundary condition: CORE bulk
33   USE sbcblk_clio     ! Surface boundary condition: CLIO bulk
34   USE sbccpl          ! Surface boundary condition: coupled interface
35   USE albedo          ! ocean & ice albedo
36
37   USE phycst          ! Define parameters for the routines
38   USE eosbn2          ! equation of state
39   USE limdyn          ! Ice dynamics
40   USE limtrp          ! Ice transport
41   USE limthd          ! Ice thermodynamics
42   USE limitd_th       ! Thermodynamics on ice thickness distribution
43   USE limitd_me       ! Mechanics on ice thickness distribution
44   USE limsbc          ! sea surface boundary condition
45   USE limdiahsb       ! Ice budget diagnostics
46   USE limwri          ! Ice outputs
47   USE limrst          ! Ice restarts
48   USE limupdate1       ! update of global variables
49   USE limupdate2       ! update of global variables
50   USE limvar          ! Ice variables switch
51
52   USE c1d             ! 1D vertical configuration
53   USE lbclnk          ! lateral boundary condition - MPP link
54   USE lib_mpp         ! MPP library
55   USE wrk_nemo        ! work arrays
56   USE timing          ! Timing
57   USE iom             ! I/O manager library
58   USE in_out_manager  ! I/O manager
59   USE prtctl          ! Print control
60   USE lib_fortran     !
61
62#if defined key_bdy 
63   USE bdyice_lim       ! unstructured open boundary data  (bdy_ice_lim routine)
64#endif
65
66   IMPLICIT NONE
67   PRIVATE
68
69   PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90
70   
71   !! * Substitutions
72#  include "domzgr_substitute.h90"
73#  include "vectopt_loop_substitute.h90"
74   !!----------------------------------------------------------------------
75   !! NEMO/OPA 4.0 , UCL NEMO Consortium (2011)
76   !! $Id$
77   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
78   !!----------------------------------------------------------------------
79CONTAINS
80
81   FUNCTION fice_cell_ave ( ptab)
82      !!--------------------------------------------------------------------------
83      !! * Compute average over categories, for grid cell (ice covered and free ocean)
84      !!--------------------------------------------------------------------------
85      REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave
86      REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab
87      INTEGER :: jl ! Dummy loop index
88     
89      fice_cell_ave (:,:) = 0.0_wp
90     
91      DO jl = 1, jpl
92         fice_cell_ave (:,:) = fice_cell_ave (:,:) &
93            &                  + a_i (:,:,jl) * ptab (:,:,jl)
94      END DO
95     
96   END FUNCTION fice_cell_ave
97   
98   FUNCTION fice_ice_ave ( ptab)
99      !!--------------------------------------------------------------------------
100      !! * Compute average over categories, for ice covered part of grid cell
101      !!--------------------------------------------------------------------------
102      REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave
103      REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab
104
105      fice_ice_ave (:,:) = 0.0_wp
106      WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)
107
108   END FUNCTION fice_ice_ave
109
110   !!======================================================================
111
112   SUBROUTINE sbc_ice_lim( kt, kblk )
113      !!---------------------------------------------------------------------
114      !!                  ***  ROUTINE sbc_ice_lim  ***
115      !!                   
116      !! ** Purpose :   update the ocean surface boundary condition via the
117      !!                Louvain la Neuve Sea Ice Model time stepping
118      !!
119      !! ** Method  :   ice model time stepping
120      !!              - call the ice dynamics routine
121      !!              - call the ice advection/diffusion routine
122      !!              - call the ice thermodynamics routine
123      !!              - call the routine that computes mass and
124      !!                heat fluxes at the ice/ocean interface
125      !!              - save the outputs
126      !!              - save the outputs for restart when necessary
127      !!
128      !! ** Action  : - time evolution of the LIM sea-ice model
129      !!              - update all sbc variables below sea-ice:
130      !!                utau, vtau, taum, wndm, qns , qsr, emp , sfx
131      !!---------------------------------------------------------------------
132      INTEGER, INTENT(in) ::   kt      ! ocean time step
133      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE)
134      !!
135      INTEGER  ::   jl      ! dummy loop index
136      REAL(wp) ::   zcoef   ! local scalar
137      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky
138      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice      ! mean albedo of ice (for coupled)
139
140      REAL(wp), POINTER, DIMENSION(:,:) :: zalb_ice_all    ! Mean albedo over all categories
141      REAL(wp), POINTER, DIMENSION(:,:) :: ztem_ice_all    ! Mean temperature over all categories
142     
143      REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_ice_all   ! Mean solar heat flux over all categories
144      REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_ice_all   ! Mean non solar heat flux over all categories
145      REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_ice_all   ! Mean latent heat flux over all categories
146      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all  ! Mean d(qns)/dT over all categories
147      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqla_ice_all  ! Mean d(qla)/dT over all categories
148      !!----------------------------------------------------------------------
149
150      !- O.M. : why do we allocate all these arrays even when MOD( kt-1, nn_fsbc ) /= 0 ?????
151
152      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim')
153
154      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
155
156#if defined key_coupled
157      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_alloc( jpi,jpj,jpl, zalb_ice)
158      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
159         &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
160#endif
161
162      IF( kt == nit000 ) THEN
163         IF(lwp) WRITE(numout,*)
164         IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition' 
165         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping'
166         !
167         CALL ice_init
168         !
169         IF( ln_nicep ) THEN      ! control print at a given point
170            jiindx = 177   ;   jjindx = 112
171            IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
172         ENDIF
173      ENDIF
174
175      !                                        !----------------------!
176      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  !
177         !                                     !----------------------!
178         !                                           !  Bulk Formulea !
179         !                                           !----------------!
180         !
181         u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point
182         v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean)
183         !
184         t_bo(:,:) = tfreez( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin]
185         !                                           ! (set to rt0 over land)
186         CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo
187
188         DO jl = 1, jpl
189            t_su(:,:,jl) = t_su(:,:,jl) +  rt0 * ( 1. - tmask(:,:,1) )
190         END DO
191
192         IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) )
193         
194#if defined key_coupled
195         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
196            !
197            ! Compute mean albedo and temperature
198            zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) ) 
199            ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) ) 
200            !
201         ENDIF
202#endif
203                                               ! Bulk formulea - provides the following fields:
204         ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2]
205         ! qsr_ice , qns_ice  : solar & non solar heat flux over ice   (T-point)         [W/m2]
206         ! qla_ice            : latent heat flux over ice              (T-point)         [W/m2]
207         ! dqns_ice, dqla_ice : non solar & latent heat sensistivity   (T-point)         [W/m2]
208         ! tprecip , sprecip  : total & solid precipitation            (T-point)         [Kg/m2/s]
209         ! fr1_i0  , fr2_i0   : 1sr & 2nd fraction of qsr penetration in ice             [%]
210         !
211         SELECT CASE( kblk )
212         CASE( 3 )                                       ! CLIO bulk formulation
213            CALL blk_ice_clio( t_su , zalb_ice_cs, zalb_ice_os,                           &
214               &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   &
215               &                      qla_ice    , dqns_ice   , dqla_ice  ,               &
216               &                      tprecip    , sprecip    ,                           &
217               &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  )
218            !         
219         CASE( 4 )                                       ! CORE bulk formulation
220            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice_cs,               &
221               &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   &
222               &                      qla_ice   , dqns_ice  , dqla_ice   ,               &
223               &                      tprecip   , sprecip   ,                            &
224               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  )
225            !
226         CASE ( 5 )
227            zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) )
228           
229            CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )
230
231            CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=tn_ice    )
232
233            ! Latent heat flux is forced to 0 in coupled :
234            !  it is included in qns (non-solar heat flux)
235            qla_ice  (:,:,:) = 0.0e0_wp
236            dqla_ice (:,:,:) = 0.0e0_wp
237            !
238         END SELECT
239
240         ! Average over all categories
241#if defined key_coupled
242         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
243
244            z_qns_ice_all  (:,:) = fice_ice_ave ( qns_ice  (:,:,:) )
245            z_qsr_ice_all  (:,:) = fice_ice_ave ( qsr_ice  (:,:,:) )
246            z_dqns_ice_all (:,:) = fice_ice_ave ( dqns_ice (:,:,:) )
247            z_qla_ice_all  (:,:) = fice_ice_ave ( qla_ice  (:,:,:) )
248            z_dqla_ice_all (:,:) = fice_ice_ave ( dqla_ice (:,:,:) )
249
250            DO jl = 1, jpl
251               dqns_ice (:,:,jl) = z_dqns_ice_all (:,:)
252               dqla_ice (:,:,jl) = z_dqla_ice_all (:,:)
253            END DO
254            !
255            IF ( ln_iceflx_ave ) THEN
256               DO jl = 1, jpl
257                  qns_ice  (:,:,jl) = z_qns_ice_all  (:,:)
258                  qsr_ice  (:,:,jl) = z_qsr_ice_all  (:,:)
259                  qla_ice  (:,:,jl) = z_qla_ice_all  (:,:)
260               END DO
261            END IF
262            !
263            IF ( ln_iceflx_linear ) THEN
264               DO jl = 1, jpl
265                  qns_ice  (:,:,jl) = z_qns_ice_all(:,:) + z_dqns_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:))
266                  qla_ice  (:,:,jl) = z_qla_ice_all(:,:) + z_dqla_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:))
267                  qsr_ice  (:,:,jl) = (1.0e0_wp-zalb_ice(:,:,jl)) / (1.0e0_wp-zalb_ice_all(:,:)) * z_qsr_ice_all(:,:)
268               END DO
269            END IF
270         END IF
271#endif
272         !                                           !----------------------!
273         !                                           ! LIM-3  time-stepping !
274         !                                           !----------------------!
275         !
276         numit = numit + nn_fsbc                     ! Ice model time step
277         !
278         !                                           ! Store previous ice values
279!!gm : remark   old_...   should becomes ...b  as tn versus tb 
280         old_a_i  (:,:,:)   = a_i  (:,:,:)     ! ice area
281         old_e_i  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy
282         old_v_i  (:,:,:)   = v_i  (:,:,:)     ! ice volume
283         old_v_s  (:,:,:)   = v_s  (:,:,:)     ! snow volume
284         old_e_s  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy
285         old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content
286         old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content
287         !
288         old_u_ice(:,:) = u_ice(:,:)
289         old_v_ice(:,:) = v_ice(:,:)
290         !                                           ! intialisation to zero    !!gm is it truly necessary ???
291         d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp
292         d_v_i_thd  (:,:,:)   = 0._wp   ;   d_v_i_trp  (:,:,:)   = 0._wp
293         d_e_i_thd  (:,:,:,:) = 0._wp   ;   d_e_i_trp  (:,:,:,:) = 0._wp
294         d_v_s_thd  (:,:,:)   = 0._wp   ;   d_v_s_trp  (:,:,:)   = 0._wp
295         d_e_s_thd  (:,:,:,:) = 0._wp   ;   d_e_s_trp  (:,:,:,:) = 0._wp
296         d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp
297         d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp
298         !
299         d_u_ice_dyn(:,:) = 0._wp
300         d_v_ice_dyn(:,:) = 0._wp
301         !
302         sfx    (:,:) = 0._wp   ;   sfx_thd  (:,:) = 0._wp
303         sfx_bri(:,:) = 0._wp   ;   sfx_mec  (:,:) = 0._wp   ;   sfx_res  (:,:) = 0._wp
304         fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp
305         fhmec  (:,:) = 0._wp   ;   
306         fmmec  (:,:) = 0._wp
307         fmmflx (:,:) = 0._wp     
308         focea2D(:,:) = 0._wp
309         fsup2D (:,:) = 0._wp
310
311         ! used in limthd.F90
312         rdvosif(:,:) = 0._wp   ! variation of ice volume at surface
313         rdvobif(:,:) = 0._wp   ! variation of ice volume at bottom
314         fdvolif(:,:) = 0._wp   ! total variation of ice volume
315         rdvonif(:,:) = 0._wp   ! lateral variation of ice volume
316         fstric (:,:) = 0._wp   ! part of solar radiation transmitted through the ice
317         ffltbif(:,:) = 0._wp   ! linked with fstric
318         qfvbq  (:,:) = 0._wp   ! linked with fstric
319         rdm_snw(:,:) = 0._wp   ! variation of snow mass per unit area
320         rdm_ice(:,:) = 0._wp   ! variation of ice mass per unit area
321         hicifp (:,:) = 0._wp   ! daily thermodynamic ice production.
322         !
323         diag_sni_gr(:,:) = 0._wp   ;   diag_lat_gr(:,:) = 0._wp
324         diag_bot_gr(:,:) = 0._wp   ;   diag_dyn_gr(:,:) = 0._wp
325         diag_bot_me(:,:) = 0._wp   ;   diag_sur_me(:,:) = 0._wp
326         diag_res_pr(:,:) = 0._wp   ;   diag_trp_vi(:,:) = 0._wp
327         ! dynamical invariants
328         delta_i(:,:) = 0._wp       ;   divu_i(:,:) = 0._wp       ;   shear_i(:,:) = 0._wp
329
330                          CALL lim_rst_opn( kt )     ! Open Ice restart file
331         !
332         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print
333         ! ----------------------------------------------
334         ! ice dynamics and transport (except in 1D case)
335         ! ----------------------------------------------
336         IF( .NOT. lk_c1d ) THEN
337                          CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics )
338                          CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion )
339                          CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting
340         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print
341                          CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting)
342                          CALL lim_var_agg( 1 ) 
343#if defined key_bdy
344                          ! bdy ice thermo
345                          CALL lim_var_glo2eqv            ! equivalent variables
346                          CALL bdy_ice_lim( kt )
347                          CALL lim_itd_me_zapsmall
348                          CALL lim_var_agg(1)
349         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print
350#endif
351                          CALL lim_update1
352         ENDIF
353!                         !- Change old values for new values
354                          old_u_ice(:,:)   = u_ice (:,:)
355                          old_v_ice(:,:)   = v_ice (:,:)
356                          old_a_i(:,:,:)   = a_i (:,:,:)
357                          old_v_s(:,:,:)   = v_s (:,:,:)
358                          old_v_i(:,:,:)   = v_i (:,:,:)
359                          old_e_s(:,:,:,:) = e_s (:,:,:,:)
360                          old_e_i(:,:,:,:) = e_i (:,:,:,:)
361                          old_oa_i(:,:,:)  = oa_i(:,:,:)
362                          old_smv_i(:,:,:) = smv_i (:,:,:)
363 
364         ! ----------------------------------------------
365         ! ice thermodynamic
366         ! ----------------------------------------------
367                          CALL lim_var_glo2eqv            ! equivalent variables
368                          CALL lim_var_agg(1)             ! aggregate ice categories
369                          ! previous lead fraction and ice volume for flux calculations
370                          pfrld(:,:)   = 1._wp - at_i(:,:)
371                          phicif(:,:)  = vt_i(:,:)
372                          !
373                          CALL lim_var_bv                 ! bulk brine volume (diag)
374                          CALL lim_thd( kt )              ! Ice thermodynamics
375                          zcoef = rdt_ice /rday           !  Ice natural aging
376                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef
377                          CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin)
378         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print
379                          CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  !
380                          CALL lim_var_agg( 1 )           ! requested by limupdate
381                          CALL lim_update2                ! Global variables update
382
383                          CALL lim_var_glo2eqv            ! equivalent variables (outputs)
384                          CALL lim_var_agg(2)             ! aggregate ice thickness categories
385         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print
386         !
387                          CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes
388         !
389         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print
390         !
391         !                                           ! Diagnostics and outputs
392         IF (ln_limdiaout) CALL lim_diahsb
393!clem # if ! defined key_iomput
394                          CALL lim_wri( 1  )              ! Ice outputs
395!clem # endif
396         IF( kt == nit000 .AND. ln_rstart )   &
397            &             CALL iom_close( numrir )        ! clem: close input ice restart file
398         !
399         IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file
400                          CALL lim_var_glo2eqv            ! ???
401         !
402         IF( ln_nicep )   CALL lim_ctl( kt )              ! alerts in case of model crash
403         !
404      ENDIF                                    ! End sea-ice time step only
405
406      !                                        !--------------------------!
407      !                                        !  at all ocean time step  !
408      !                                        !--------------------------!
409      !                                               
410      !                                              ! Update surface ocean stresses (only in ice-dynamic case)
411      !                                                   ! otherwise the atm.-ocean stresses are used everywhere
412      IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents
413     
414!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
415      !
416      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
417
418#if defined key_coupled
419      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice)
420      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
421         &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
422#endif
423      !
424      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim')
425      !
426   END SUBROUTINE sbc_ice_lim
427
428
429   SUBROUTINE lim_ctl( kt )
430      !!-----------------------------------------------------------------------
431      !!                   ***  ROUTINE lim_ctl ***
432      !!                 
433      !! ** Purpose :   Alerts in case of model crash
434      !!-------------------------------------------------------------------
435      INTEGER, INTENT(in) ::   kt      ! ocean time step
436      INTEGER  ::   ji, jj, jk,  jl   ! dummy loop indices
437      INTEGER  ::   inb_altests       ! number of alert tests (max 20)
438      INTEGER  ::   ialert_id         ! number of the current alert
439      REAL(wp) ::   ztmelts           ! ice layer melting point
440      CHARACTER (len=30), DIMENSION(20)      ::   cl_alname   ! name of alert
441      INTEGER           , DIMENSION(20)      ::   inb_alp     ! number of alerts positive
442      !!-------------------------------------------------------------------
443
444      inb_altests = 10
445      inb_alp(:)  =  0
446
447      ! Alert if incompatible volume and concentration
448      ialert_id = 2 ! reference number of this alert
449      cl_alname(ialert_id) = ' Incompat vol and con         '    ! name of the alert
450
451      DO jl = 1, jpl
452         DO jj = 1, jpj
453            DO ji = 1, jpi
454               IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN
455                  !WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration '
456                  !WRITE(numout,*) ' at_i     ', at_i(ji,jj)
457                  !WRITE(numout,*) ' Point - category', ji, jj, jl
458                  !WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl)
459                  !WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl)
460                  !WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
461                  !WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
462                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1
463               ENDIF
464            END DO
465         END DO
466      END DO
467
468      ! Alerte if very thick ice
469      ialert_id = 3 ! reference number of this alert
470      cl_alname(ialert_id) = ' Very thick ice               ' ! name of the alert
471      jl = jpl 
472      DO jj = 1, jpj
473         DO ji = 1, jpi
474            IF(   ht_i(ji,jj,jl)  >  50._wp   ) THEN
475               !CALL lim_prt_state( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' )
476               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
477            ENDIF
478         END DO
479      END DO
480
481      ! Alert if very fast ice
482      ialert_id = 4 ! reference number of this alert
483      cl_alname(ialert_id) = ' Very fast ice               ' ! name of the alert
484      DO jj = 1, jpj
485         DO ji = 1, jpi
486            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5  .AND.  &
487               &  at_i(ji,jj) > 0._wp   ) THEN
488               !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 4 :   Very fast ice ' )
489               !WRITE(numout,*) ' ice strength             : ', strength(ji,jj)
490               !WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj)
491               !WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj)
492               !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj)
493               !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj)
494               !WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj)
495               !WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj)
496               !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj)
497               !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj)
498               !WRITE(numout,*)
499               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
500            ENDIF
501         END DO
502      END DO
503
504      ! Alert if there is ice on continents
505      ialert_id = 6 ! reference number of this alert
506      cl_alname(ialert_id) = ' Ice on continents           ' ! name of the alert
507      DO jj = 1, jpj
508         DO ji = 1, jpi
509            IF(   tms(ji,jj) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN 
510               !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 6 :   Ice on continents ' )
511               !WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)
512               !WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
513               !WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
514               !WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj)
515               !WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj)
516               !WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1)
517               !WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj)
518               !WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj)
519               !
520               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
521            ENDIF
522         END DO
523      END DO
524
525!
526!     ! Alert if very fresh ice
527      ialert_id = 7 ! reference number of this alert
528      cl_alname(ialert_id) = ' Very fresh ice               ' ! name of the alert
529      DO jl = 1, jpl
530         DO jj = 1, jpj
531            DO ji = 1, jpi
532               IF( sm_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN
533!                 CALL lim_prt_state(kt,ji,jj,1, ' ALERTE 7 :   Very fresh ice ' )
534!                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
535!                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
536!                 WRITE(numout,*) ' s_i_newice           : ', s_i_newice(ji,jj,1:jpl)
537!                 WRITE(numout,*)
538                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1
539               ENDIF
540            END DO
541         END DO
542      END DO
543!
544
545!     ! Alert if too old ice
546      ialert_id = 9 ! reference number of this alert
547      cl_alname(ialert_id) = ' Very old   ice               ' ! name of the alert
548      DO jl = 1, jpl
549         DO jj = 1, jpj
550            DO ji = 1, jpi
551               IF ( ( ( ABS( o_i(ji,jj,jl) ) > rdt_ice ) .OR. &
552                      ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. &
553                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN
554                  !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ')
555                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1
556               ENDIF
557            END DO
558         END DO
559      END DO
560 
561      ! Alert on salt flux
562      ialert_id = 5 ! reference number of this alert
563      cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert
564      DO jj = 1, jpj
565         DO ji = 1, jpi
566            IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth
567               !CALL lim_prt_state( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' )
568               !DO jl = 1, jpl
569                  !WRITE(numout,*) ' Category no: ', jl
570                  !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)   
571                  !WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)
572                  !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)   
573                  !WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)
574                  !WRITE(numout,*) ' '
575               !END DO
576               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
577            ENDIF
578         END DO
579      END DO
580
581      ! Alert if qns very big
582      ialert_id = 8 ! reference number of this alert
583      cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert
584      DO jj = 1, jpj
585         DO ji = 1, jpi
586            IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN
587               !
588               !WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux'
589               !WRITE(numout,*) ' ji, jj    : ', ji, jj
590               !WRITE(numout,*) ' qns       : ', qns(ji,jj)
591               !WRITE(numout,*) ' sst       : ', sst_m(ji,jj)
592               !WRITE(numout,*) ' sss       : ', sss_m(ji,jj)
593               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj)
594               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj)
595               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice
596               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice
597               !WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj)
598               !WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj)
599               !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice
600               !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice
601               !WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)
602               !WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj)
603               !WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj)
604               !WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)
605               !WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj)
606               !
607               !CALL lim_prt_state( kt, ji, jj, 2, '   ')
608               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
609               !
610            ENDIF
611         END DO
612      END DO
613      !+++++
614 
615      ! Alert if very warm ice
616      ialert_id = 10 ! reference number of this alert
617      cl_alname(ialert_id) = ' Very warm ice                ' ! name of the alert
618      inb_alp(ialert_id) = 0
619      DO jl = 1, jpl
620         DO jk = 1, nlay_i
621            DO jj = 1, jpj
622               DO ji = 1, jpi
623                  ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rtt
624                  IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   &
625                     &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN
626                     !WRITE(numout,*) ' ALERTE 10 :   Very warm ice'
627                     !WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl
628                     !WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl)
629                     !WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
630                     !WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl)
631                     !WRITE(numout,*) ' ztmelts : ', ztmelts
632                     inb_alp(ialert_id) = inb_alp(ialert_id) + 1
633                  ENDIF
634               END DO
635            END DO
636         END DO
637      END DO
638
639      ! sum of the alerts on all processors
640      IF( lk_mpp ) THEN
641         DO ialert_id = 1, inb_altests
642            CALL mpp_sum(inb_alp(ialert_id))
643         END DO
644      ENDIF
645
646      ! print alerts
647      IF( lwp ) THEN
648         ialert_id = 1                                 ! reference number of this alert
649         cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert
650         WRITE(numout,*) ' time step ',kt
651         WRITE(numout,*) ' All alerts at the end of ice model '
652         DO ialert_id = 1, inb_altests
653            WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! '
654         END DO
655      ENDIF
656     !
657   END SUBROUTINE lim_ctl
658 
659   
660   SUBROUTINE lim_prt_state( kt, ki, kj, kn, cd1 )
661      !!-----------------------------------------------------------------------
662      !!                   ***  ROUTINE lim_prt_state ***
663      !!                 
664      !! ** Purpose :   Writes global ice state on the (i,j) point
665      !!                in ocean.ouput
666      !!                3 possibilities exist
667      !!                n = 1/-1 -> simple ice state (plus Mechanical Check if -1)
668      !!                n = 2    -> exhaustive state
669      !!                n = 3    -> ice/ocean salt fluxes
670      !!
671      !! ** input   :   point coordinates (i,j)
672      !!                n : number of the option
673      !!-------------------------------------------------------------------
674      INTEGER         , INTENT(in) ::   kt      ! ocean time step
675      INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices
676      CHARACTER(len=*), INTENT(in) ::   cd1           !
677      !!
678      INTEGER :: jl, ji, jj
679      !!-------------------------------------------------------------------
680
681      DO ji = mi0(ki), mi1(ki)
682         DO jj = mj0(kj), mj1(kj)
683
684            WRITE(numout,*) ' time step ',kt,' ',cd1             ! print title
685
686            !----------------
687            !  Simple state
688            !----------------
689           
690            IF ( kn == 1 .OR. kn == -1 ) THEN
691               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj
692               WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
693               WRITE(numout,*) ' Simple state '
694               WRITE(numout,*) ' masks s,u,v   : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)
695               WRITE(numout,*) ' lat - long    : ', gphit(ji,jj), glamt(ji,jj)
696               WRITE(numout,*) ' Time step     : ', numit
697               WRITE(numout,*) ' - Ice drift   '
698               WRITE(numout,*) '   ~~~~~~~~~~~ '
699               WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj)
700               WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj)
701               WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1)
702               WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj)
703               WRITE(numout,*) ' strength      : ', strength(ji,jj)
704               WRITE(numout,*)
705               WRITE(numout,*) ' - Cell values '
706               WRITE(numout,*) '   ~~~~~~~~~~~ '
707               WRITE(numout,*) ' cell area     : ', area(ji,jj)
708               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)       
709               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)       
710               WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)       
711               DO jl = 1, jpl
712                  WRITE(numout,*) ' - Category (', jl,')'
713                  WRITE(numout,*) ' a_i           : ', a_i(ji,jj,jl)
714                  WRITE(numout,*) ' ht_i          : ', ht_i(ji,jj,jl)
715                  WRITE(numout,*) ' ht_s          : ', ht_s(ji,jj,jl)
716                  WRITE(numout,*) ' v_i           : ', v_i(ji,jj,jl)
717                  WRITE(numout,*) ' v_s           : ', v_s(ji,jj,jl)
718                  WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl)/1.0e9
719                  WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl)/1.0e9
720                  WRITE(numout,*) ' t_su          : ', t_su(ji,jj,jl)
721                  WRITE(numout,*) ' t_snow        : ', t_s(ji,jj,1,jl)
722                  WRITE(numout,*) ' t_i           : ', t_i(ji,jj,1:nlay_i,jl)
723                  WRITE(numout,*) ' sm_i          : ', sm_i(ji,jj,jl)
724                  WRITE(numout,*) ' smv_i         : ', smv_i(ji,jj,jl)
725                  WRITE(numout,*)
726               END DO
727            ENDIF
728            IF( kn == -1 ) THEN
729               WRITE(numout,*) ' Mechanical Check ************** '
730               WRITE(numout,*) ' Check what means ice divergence '
731               WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj)
732               WRITE(numout,*) ' Total lead fraction     ', ato_i(ji,jj)
733               WRITE(numout,*) ' Sum of both             ', ato_i(ji,jj) + at_i(ji,jj)
734               WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ji,jj) + at_i(ji,jj) - 1.00
735            ENDIF
736           
737
738            !--------------------
739            !  Exhaustive state
740            !--------------------
741           
742            IF ( kn .EQ. 2 ) THEN
743               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj
744               WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
745               WRITE(numout,*) ' Exhaustive state '
746               WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
747               WRITE(numout,*) ' Time step ', numit
748               WRITE(numout,*) 
749               WRITE(numout,*) ' - Cell values '
750               WRITE(numout,*) '   ~~~~~~~~~~~ '
751               WRITE(numout,*) ' cell area     : ', area(ji,jj)
752               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)       
753               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)       
754               WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)       
755               WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj)
756               WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj)
757               WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1)
758               WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj)
759               WRITE(numout,*) ' strength      : ', strength(ji,jj)
760               WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ji,jj)
761               WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ji,jj)  , ' old_v_ice     : ', old_v_ice(ji,jj) 
762               WRITE(numout,*)
763               
764               DO jl = 1, jpl
765                  WRITE(numout,*) ' - Category (',jl,')'
766                  WRITE(numout,*) '   ~~~~~~~~         ' 
767                  WRITE(numout,*) ' ht_i       : ', ht_i(ji,jj,jl)             , ' ht_s       : ', ht_s(ji,jj,jl)
768                  WRITE(numout,*) ' t_i        : ', t_i(ji,jj,1:nlay_i,jl)
769                  WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1,jl)
770                  WRITE(numout,*) ' sm_i       : ', sm_i(ji,jj,jl)             , ' o_i        : ', o_i(ji,jj,jl)
771                  WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' old_a_i    : ', old_a_i(ji,jj,jl)   
772                  WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl) 
773                  WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' old_v_i    : ', old_v_i(ji,jj,jl)   
774                  WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl) 
775                  WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' old_v_s    : ', old_v_s(ji,jj,jl) 
776                  WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ji,jj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ji,jj,jl)
777                  WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ji,jj,1,jl)/1.0e9 
778                  WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ji,jj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ji,jj,1,jl)/1.0e9
779                  WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ji,jj,2,jl)/1.0e9 
780                  WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ji,jj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ji,jj,2,jl)/1.0e9
781                  WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' old_e_snow : ', old_e_s(ji,jj,1,jl) 
782                  WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ji,jj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ji,jj,1,jl)
783                  WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' old_smv_i  : ', old_smv_i(ji,jj,jl)   
784                  WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl) 
785                  WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' old_oa_i   : ', old_oa_i(ji,jj,jl)
786                  WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl)
787               END DO !jl
788               
789               WRITE(numout,*)
790               WRITE(numout,*) ' - Heat / FW fluxes '
791               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
792               WRITE(numout,*) ' emp        : ', emp      (ji,jj)
793               WRITE(numout,*) ' sfx        : ', sfx      (ji,jj)
794               WRITE(numout,*) ' sfx_thd    : ', sfx_thd(ji,jj)
795               WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ji,jj)
796               WRITE(numout,*) ' sfx_mec    : ', sfx_mec  (ji,jj)
797               WRITE(numout,*) ' sfx_res    : ', sfx_res(ji,jj)
798               WRITE(numout,*) ' fmmec      : ', fmmec    (ji,jj)
799               WRITE(numout,*) ' fhmec      : ', fhmec    (ji,jj)
800               WRITE(numout,*) ' fhbri      : ', fhbri    (ji,jj)
801               WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ji,jj)
802               WRITE(numout,*) 
803               WRITE(numout,*) ' sst        : ', sst_m(ji,jj) 
804               WRITE(numout,*) ' sss        : ', sss_m(ji,jj) 
805               WRITE(numout,*) 
806               WRITE(numout,*) ' - Stresses '
807               WRITE(numout,*) '   ~~~~~~~~ '
808               WRITE(numout,*) ' utau_ice   : ', utau_ice(ji,jj) 
809               WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ji,jj)
810               WRITE(numout,*) ' utau       : ', utau    (ji,jj) 
811               WRITE(numout,*) ' vtau       : ', vtau    (ji,jj)
812               WRITE(numout,*) ' oc. vel. u : ', u_oce   (ji,jj)
813               WRITE(numout,*) ' oc. vel. v : ', v_oce   (ji,jj)
814            ENDIF
815           
816            !---------------------
817            ! Salt / heat fluxes
818            !---------------------
819           
820            IF ( kn .EQ. 3 ) THEN
821               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj
822               WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
823               WRITE(numout,*) ' - Salt / Heat Fluxes '
824               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
825               WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
826               WRITE(numout,*) ' Time step ', numit
827               WRITE(numout,*)
828               WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
829               WRITE(numout,*) ' qsr       : ', qsr(ji,jj)
830               WRITE(numout,*) ' qns       : ', qns(ji,jj)
831               WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)
832               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) * r1_rdtice
833               WRITE(numout,*) ' qldif     : ', qldif(ji,jj) * r1_rdtice
834               WRITE(numout,*)
835               WRITE(numout,*) ' - Salt fluxes at bottom interface ***'
836               WRITE(numout,*) ' emp       : ', emp    (ji,jj)
837               WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj)
838               WRITE(numout,*) ' sfx       : ', sfx    (ji,jj)
839               WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj)
840               WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ji,jj)
841               WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
842               WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)
843               WRITE(numout,*)
844               WRITE(numout,*) ' - Momentum fluxes '
845               WRITE(numout,*) ' utau      : ', utau(ji,jj) 
846               WRITE(numout,*) ' vtau      : ', vtau(ji,jj)
847            ENDIF
848            WRITE(numout,*) ' '
849            !
850         END DO
851      END DO
852
853   END SUBROUTINE lim_prt_state
854
855#else
856   !!----------------------------------------------------------------------
857   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
858   !!----------------------------------------------------------------------
859CONTAINS
860   SUBROUTINE sbc_ice_lim ( kt, kblk )     ! Dummy routine
861      WRITE(*,*) 'sbc_ice_lim: You should not have seen this print! error?', kt, kblk
862   END SUBROUTINE sbc_ice_lim
863#endif
864
865   !!======================================================================
866END MODULE sbcice_lim
Note: See TracBrowser for help on using the repository browser.