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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 35.3 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   !!            4.0  ! 2011-01  (A Porter)  dynamical allocation
13   !!----------------------------------------------------------------------
14#if defined key_lim3
15   !!----------------------------------------------------------------------
16   !!   'key_lim3' :                                  LIM 3.0 sea-ice model
17   !!----------------------------------------------------------------------
18   !!   sbc_ice_lim  : sea-ice model time-stepping and update ocean sbc over ice-covered area
19   !!   lim_ctl       : alerts in case of ice model crash
20   !!   lim_prt_state : ice control print at a given grid point
21   !!----------------------------------------------------------------------
22   USE oce             ! ocean dynamics and tracers
23   USE dom_oce         ! ocean space and time domain
24   USE par_ice         ! sea-ice parameters
25   USE ice             ! LIM-3: ice variables
26   USE iceini          ! LIM-3: ice initialisation
27   USE dom_ice         ! LIM-3: ice domain
28
29   USE sbc_oce         ! Surface boundary condition: ocean fields
30   USE sbc_ice         ! Surface boundary condition: ice   fields
31   USE sbcblk_core     ! Surface boundary condition: CORE bulk
32   USE sbcblk_clio     ! Surface boundary condition: CLIO bulk
33   USE albedo          ! ocean & ice albedo
34
35   USE phycst          ! Define parameters for the routines
36   USE eosbn2          ! equation of state
37   USE limdyn          ! Ice dynamics
38   USE limtrp          ! Ice transport
39   USE limthd          ! Ice thermodynamics
40   USE limitd_th       ! Thermodynamics on ice thickness distribution
41   USE limitd_me       ! Mechanics on ice thickness distribution
42   USE limsbc          ! sea surface boundary condition
43   USE limdia          ! Ice diagnostics
44   USE limwri          ! Ice outputs
45   USE limrst          ! Ice restarts
46   USE limupdate       ! update of global variables
47   USE limvar          ! Ice variables switch
48
49   USE c1d             ! 1D vertical configuration
50   USE lbclnk          ! lateral boundary condition - MPP link
51   USE lib_mpp         ! MPP library
52   USE iom             ! I/O manager library
53   USE in_out_manager  ! I/O manager
54   USE prtctl          ! Print control
55
56   IMPLICIT NONE
57   PRIVATE
58
59   PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90
60
61   !! * Control permutation of array indices
62#  include "oce_ftrans.h90"
63#  include "dom_oce_ftrans.h90"
64#  include "sbc_oce_ftrans.h90"
65   
66   !! * Substitutions
67#  include "domzgr_substitute.h90"
68#  include "vectopt_loop_substitute.h90"
69   !!----------------------------------------------------------------------
70   !! NEMO/OPA 4.0 , UCL NEMO Consortium (2011)
71   !! $Id$
72   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
73   !!----------------------------------------------------------------------
74CONTAINS
75
76   SUBROUTINE sbc_ice_lim( kt, kblk )
77      !!---------------------------------------------------------------------
78      !!                  ***  ROUTINE sbc_ice_lim  ***
79      !!                   
80      !! ** Purpose :   update the ocean surface boundary condition via the
81      !!                Louvain la Neuve Sea Ice Model time stepping
82      !!
83      !! ** Method  :   ice model time stepping
84      !!              - call the ice dynamics routine
85      !!              - call the ice advection/diffusion routine
86      !!              - call the ice thermodynamics routine
87      !!              - call the routine that computes mass and
88      !!                heat fluxes at the ice/ocean interface
89      !!              - save the outputs
90      !!              - save the outputs for restart when necessary
91      !!
92      !! ** Action  : - time evolution of the LIM sea-ice model
93      !!              - update all sbc variables below sea-ice:
94      !!                utau, vtau, taum, wndm, qns , qsr, emp , emps
95      !!---------------------------------------------------------------------
96      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
97      USE wrk_nemo, ONLY:   zalb_ice_os => wrk_3d_1 ! albedo of the ice under overcast sky
98      USE wrk_nemo, ONLY:   zalb_ice_cs => wrk_3d_2 ! albedo of ice under clear sky
99      !!
100      INTEGER, INTENT(in) ::   kt      ! ocean time step
101      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE)
102      !!
103      INTEGER  ::   jl      ! dummy loop index
104      REAL(wp) ::   zcoef   ! local scalar
105      !!----------------------------------------------------------------------
106
107      IF( wrk_in_use(3, 1,2) ) THEN
108         CALL ctl_stop( 'sbc_ice_lim: requested workspace arrays are unavailable' )   ;   RETURN
109      ENDIF
110
111      IF( kt == nit000 ) THEN
112         IF(lwp) WRITE(numout,*)
113         IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition' 
114         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping'
115         !
116         CALL ice_init
117         !
118         IF( ln_nicep ) THEN      ! control print at a given point
119            jiindx = 44   ;   jjindx = 140
120            WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
121         ENDIF
122      ENDIF
123
124      !                                        !----------------------!
125      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  !
126         !                                     !----------------------!
127         !                                           !  Bulk Formulea !
128         !                                           !----------------!
129         !
130         u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point
131         v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean)
132         !
133         t_bo(:,:) = tfreez( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin]
134         !                                           ! (set to rt0 over land)
135         CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo
136
137         DO jl = 1, jpl
138            t_su(:,:,jl) = t_su(:,:,jl) +  rt0 * ( 1. - tmask(:,:,1) )
139         END DO
140                                                     ! Bulk formulea - provides the following fields:
141         ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2]
142         ! qsr_ice , qns_ice  : solar & non solar heat flux over ice   (T-point)         [W/m2]
143         ! qla_ice            : latent heat flux over ice              (T-point)         [W/m2]
144         ! dqns_ice, dqla_ice : non solar & latent heat sensistivity   (T-point)         [W/m2]
145         ! tprecip , sprecip  : total & solid precipitation            (T-point)         [Kg/m2/s]
146         ! fr1_i0  , fr2_i0   : 1sr & 2nd fraction of qsr penetration in ice             [%]
147         !
148         SELECT CASE( kblk )
149         CASE( 3 )                                       ! CLIO bulk formulation
150            CALL blk_ice_clio( t_su , zalb_ice_cs, zalb_ice_os,                           &
151               &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   &
152               &                      qla_ice    , dqns_ice   , dqla_ice  ,               &
153               &                      tprecip    , sprecip    ,                           &
154               &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  )
155            !         
156         CASE( 4 )                                       ! CORE bulk formulation
157            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice_cs,               &
158               &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   &
159               &                      qla_ice   , dqns_ice  , dqla_ice   ,               &
160               &                      tprecip   , sprecip   ,                            &
161               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  )
162         END SELECT
163
164         !                                           !----------------------!
165         !                                           ! LIM-3  time-stepping !
166         !                                           !----------------------!
167         !
168         numit = numit + nn_fsbc                     ! Ice model time step
169         !
170         !                                           ! Store previous ice values
171!!gm : remark   old_...   should becomes ...b  as tn versus tb 
172         old_a_i  (:,:,:)   = a_i  (:,:,:)     ! ice area
173         old_e_i  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy
174         old_v_i  (:,:,:)   = v_i  (:,:,:)     ! ice volume
175         old_v_s  (:,:,:)   = v_s  (:,:,:)     ! snow volume
176         old_e_s  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy
177         old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content
178         old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content
179
180         !                                           ! intialisation to zero    !!gm is it truly necessary ???
181         d_a_i_thd  (:,:,:)   = 0.e0   ;   d_a_i_trp  (:,:,:)   = 0.e0
182         d_v_i_thd  (:,:,:)   = 0.e0   ;   d_v_i_trp  (:,:,:)   = 0.e0
183         d_e_i_thd  (:,:,:,:) = 0.e0   ;   d_e_i_trp  (:,:,:,:) = 0.e0
184         d_v_s_thd  (:,:,:)   = 0.e0   ;   d_v_s_trp  (:,:,:)   = 0.e0
185         d_e_s_thd  (:,:,:,:) = 0.e0   ;   d_e_s_trp  (:,:,:,:) = 0.e0
186         d_smv_i_thd(:,:,:)   = 0.e0   ;   d_smv_i_trp(:,:,:)   = 0.e0
187         d_oa_i_thd (:,:,:)   = 0.e0   ;   d_oa_i_trp (:,:,:)   = 0.e0
188         !
189         fseqv    (:,:) = 0.e0
190         fsbri    (:,:) = 0.e0     ;   fsalt_res(:,:) = 0.e0
191         fsalt_rpo(:,:) = 0.e0
192         fhmec    (:,:) = 0.e0     ;   fhbri    (:,:) = 0.e0
193         fmmec    (:,:) = 0.e0     ;   fheat_res(:,:) = 0.e0
194         fheat_rpo(:,:) = 0.e0     ;   focea2D  (:,:) = 0.e0
195         fsup2D   (:,:) = 0.e0
196         !
197         diag_sni_gr(:,:) = 0.e0   ;   diag_lat_gr(:,:) = 0.e0
198         diag_bot_gr(:,:) = 0.e0   ;   diag_dyn_gr(:,:) = 0.e0
199         diag_bot_me(:,:) = 0.e0   ;   diag_sur_me(:,:) = 0.e0
200         ! dynamical invariants
201         delta_i(:,:) = 0.e0       ;   divu_i(:,:) = 0.e0       ;   shear_i(:,:) = 0.e0
202
203                          CALL lim_rst_opn( kt )     ! Open Ice restart file
204         !
205         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print
206         !
207         IF( .NOT. lk_c1d ) THEN
208                                                     ! Ice dynamics & transport (not in 1D case)
209                          CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics )
210                          CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion )
211                          CALL lim_var_agg(1)             ! aggregate categories, requested
212                          CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting
213         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print
214                          CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting)
215         ENDIF
216         !                                           ! Ice thermodynamics
217                          CALL lim_var_glo2eqv            ! equivalent variables
218                          CALL lim_var_agg(1)             ! aggregate ice categories
219                          CALL lim_var_bv                 ! bulk brine volume (diag)
220                          CALL lim_thd( kt )              ! Ice thermodynamics
221                          zcoef = rdt_ice / 86400.e0      !  Ice natural aging
222                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef
223                          CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin)
224         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print
225                          CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  !
226         !
227         !                                           ! Global variables update
228                          CALL lim_var_agg( 1 )           ! requested by limupdate
229                          CALL lim_update                 ! Global variables update
230                          CALL lim_var_glo2eqv            ! equivalent variables (outputs)
231                          CALL lim_var_agg(2)             ! aggregate ice thickness categories
232         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 2, ' - Final state - ' )   ! control print
233         !
234                          CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes
235         !
236         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print
237         !
238         !                                           ! Diagnostics and outputs
239         IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   &
240            &             CALL lim_dia 
241                          CALL lim_wri( 1  )              ! Ice outputs
242         IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file
243                          CALL lim_var_glo2eqv            ! ???
244         !
245         IF( ln_nicep )   CALL lim_ctl               ! alerts in case of model crash
246         !
247      ENDIF                                    ! End sea-ice time step only
248
249      !                                        !--------------------------!
250      !                                        !  at all ocean time step  !
251      !                                        !--------------------------!
252      !                                               
253      !                                              ! Update surface ocean stresses (only in ice-dynamic case)
254      !                                                   ! otherwise the atm.-ocean stresses are used everywhere
255      IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents
256     
257!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
258      !
259      IF( wrk_not_released(3, 1,2) )   CALL ctl_stop( 'sbc_ice_lim: failed to release workspace arrays' )
260      !
261   END SUBROUTINE sbc_ice_lim
262
263
264   SUBROUTINE lim_ctl
265      !!-----------------------------------------------------------------------
266      !!                   ***  ROUTINE lim_ctl ***
267      !!                 
268      !! ** Purpose :   Alerts in case of model crash
269      !!-------------------------------------------------------------------
270      INTEGER  ::   ji, jj, jk,  jl   ! dummy loop indices
271      INTEGER  ::   inb_altests       ! number of alert tests (max 20)
272      INTEGER  ::   ialert_id         ! number of the current alert
273      REAL(wp) ::   ztmelts           ! ice layer melting point
274      CHARACTER (len=30), DIMENSION(20)      ::   cl_alname   ! name of alert
275      INTEGER           , DIMENSION(20)      ::   inb_alp     ! number of alerts positive
276      !!-------------------------------------------------------------------
277
278      inb_altests = 10
279      inb_alp(:)  = 0
280
281      ! Alert if incompatible volume and concentration
282      ialert_id = 2 ! reference number of this alert
283      cl_alname(ialert_id) = ' Incompat vol and con         '    ! name of the alert
284
285      DO jl = 1, jpl
286         DO jj = 1, jpj
287            DO ji = 1, jpi
288               IF(  v_i(ji,jj,jl) /= 0.e0   .AND.   a_i(ji,jj,jl) == 0.e0   ) THEN
289                  WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration '
290                  WRITE(numout,*) ' at_i     ', at_i(ji,jj)
291                  WRITE(numout,*) ' Point - category', ji, jj, jl
292                  WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl)
293                  WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl)
294                  WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
295                  WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
296                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1
297               ENDIF
298            END DO
299         END DO
300      END DO
301
302      ! Alerte if very thick ice
303      ialert_id = 3 ! reference number of this alert
304      cl_alname(ialert_id) = ' Very thick ice               ' ! name of the alert
305      jl = jpl 
306      DO jj = 1, jpj
307         DO ji = 1, jpi
308            IF(   ht_i(ji,jj,jl) .GT. 50.0   ) THEN
309               CALL lim_prt_state( ji, jj, 2, ' ALERTE 3 :   Very thick ice ' )
310               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
311            ENDIF
312         END DO
313      END DO
314
315      ! Alert if very fast ice
316      ialert_id = 4 ! reference number of this alert
317      cl_alname(ialert_id) = ' Very fast ice               ' ! name of the alert
318      DO jj = 1, jpj
319         DO ji = 1, jpi
320            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) .GT. 0.5  .AND.  &
321               &  at_i(ji,jj) .GT. 0.e0   ) THEN
322               CALL lim_prt_state( ji, jj, 1, ' ALERTE 4 :   Very fast ice ' )
323               WRITE(numout,*) ' ice strength             : ', strength(ji,jj)
324               WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj) 
325               WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj)
326               WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj) 
327               WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj)
328               WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj)
329               WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj)
330               WRITE(numout,*) ' sst                      : ', sst_m(ji,jj)
331               WRITE(numout,*) ' sss                      : ', sss_m(ji,jj)
332               WRITE(numout,*) 
333               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
334            ENDIF
335         END DO
336      END DO
337
338      ! Alert if there is ice on continents
339      ialert_id = 6 ! reference number of this alert
340      cl_alname(ialert_id) = ' Ice on continents           ' ! name of the alert
341      DO jj = 1, jpj
342         DO ji = 1, jpi
343            IF(   tms(ji,jj) .LE. 0.0   .AND.   at_i(ji,jj) .GT. 0.e0   ) THEN
344               CALL lim_prt_state( ji, jj, 1, ' ALERTE 6 :   Ice on continents ' )
345               WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 
346               WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
347               WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
348               WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj)
349               WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj)
350               WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1)
351               WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj)
352               WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj)
353               !
354               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
355            ENDIF
356         END DO
357      END DO
358
359!
360!     ! Alert if very fresh ice
361      ialert_id = 7 ! reference number of this alert
362      cl_alname(ialert_id) = ' Very fresh ice               ' ! name of the alert
363      DO jl = 1, jpl
364         DO jj = 1, jpj
365            DO ji = 1, jpi
366!!gm  test twice sm_i ...  ????  bug?
367               IF( ( ( ABS( sm_i(ji,jj,jl) ) .LT. 0.50) .OR. &
368                     ( ABS( sm_i(ji,jj,jl) ) .LT. 0.50) ) .AND. &
369                             ( a_i(ji,jj,jl) .GT. 0.e0 ) ) THEN
370!                 CALL lim_prt_state(ji,jj,1, ' ALERTE 7 :   Very fresh ice ' )
371!                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
372!                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
373!                 WRITE(numout,*) ' s_i_newice           : ', s_i_newice(ji,jj,1:jpl)
374!                 WRITE(numout,*)
375                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1
376               ENDIF
377            END DO
378         END DO
379      END DO
380!
381
382!     ! Alert if too old ice
383      ialert_id = 9 ! reference number of this alert
384      cl_alname(ialert_id) = ' Very old   ice               ' ! name of the alert
385      DO jl = 1, jpl
386         DO jj = 1, jpj
387            DO ji = 1, jpi
388               IF ( ( ( ABS( o_i(ji,jj,jl) ) .GT. rdt_ice ) .OR. &
389                      ( ABS( o_i(ji,jj,jl) ) .LT. 0.00) ) .AND. &
390                             ( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN
391                  CALL lim_prt_state( ji, jj, 1, ' ALERTE 9 :   Wrong ice age ')
392                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1
393               ENDIF
394            END DO
395         END DO
396      END DO
397 
398      ! Alert on salt flux
399      ialert_id = 5 ! reference number of this alert
400      cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert
401      DO jj = 1, jpj
402         DO ji = 1, jpi
403            IF( ABS( emps(ji,jj) ) .GT. 1.0e-2 ) THEN
404               CALL lim_prt_state( ji, jj, 3, ' ALERTE 5 :   High salt flux ' )
405               DO jl = 1, jpl
406                  WRITE(numout,*) ' Category no: ', jl
407                  WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)   
408                  WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl) 
409                  WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)   
410                  WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl) 
411                  WRITE(numout,*) ' '
412               END DO
413               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
414            ENDIF
415         END DO
416      END DO
417
418      ! Alert if qns very big
419      ialert_id = 8 ! reference number of this alert
420      cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert
421      DO jj = 1, jpj
422         DO ji = 1, jpi
423            IF(   ABS( qns(ji,jj) ) .GT. 1500.0   .AND.  ( at_i(ji,jj) .GT. 0.0 ) )  THEN
424               !
425               WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux'
426               WRITE(numout,*) ' ji, jj    : ', ji, jj
427               WRITE(numout,*) ' qns       : ', qns(ji,jj)
428               WRITE(numout,*) ' sst       : ', sst_m(ji,jj)
429               WRITE(numout,*) ' sss       : ', sss_m(ji,jj)
430               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj)
431               WRITE(numout,*) ' qldif     : ', qldif(ji,jj)
432               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice
433               WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice
434               WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj)
435               WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj)
436               WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice
437               WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice
438               WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj) 
439               WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj) 
440               WRITE(numout,*) ' fheat_rpo : ', fheat_rpo(ji,jj) 
441               WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 
442               WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj) 
443               !
444               CALL lim_prt_state( ji, jj, 2, '   ')
445               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
446               !
447            ENDIF
448         END DO
449      END DO
450      !+++++
451 
452      ! Alert if very warm ice
453      ialert_id = 10 ! reference number of this alert
454      cl_alname(ialert_id) = ' Very warm ice                ' ! name of the alert
455      inb_alp(ialert_id) = 0
456      DO jl = 1, jpl
457         DO jk = 1, nlay_i
458            DO jj = 1, jpj
459               DO ji = 1, jpi
460                  ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rtt
461                  IF( t_i(ji,jj,jk,jl) .GE. ztmelts  .AND.  v_i(ji,jj,jl) .GT. 1.e-6   &
462                     &                               .AND.  a_i(ji,jj,jl) .GT. 0.e0    ) THEN
463                     WRITE(numout,*) ' ALERTE 10 :   Very warm ice'
464                     WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl
465                     WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl)
466                     WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
467                     WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl)
468                     WRITE(numout,*) ' ztmelts : ', ztmelts
469                     inb_alp(ialert_id) = inb_alp(ialert_id) + 1
470                  ENDIF
471               END DO
472            END DO
473         END DO
474      END DO
475
476      ialert_id = 1                                 ! reference number of this alert
477      cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert
478      WRITE(numout,*)
479      WRITE(numout,*) ' All alerts at the end of ice model '
480      DO ialert_id = 1, inb_altests
481         WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! '
482      END DO
483      !
484   END SUBROUTINE lim_ctl
485 
486   
487   SUBROUTINE lim_prt_state( ki, kj, kn, cd1 )
488      !!-----------------------------------------------------------------------
489      !!                   ***  ROUTINE lim_prt_state ***
490      !!                 
491      !! ** Purpose :   Writes global ice state on the (i,j) point
492      !!                in ocean.ouput
493      !!                3 possibilities exist
494      !!                n = 1/-1 -> simple ice state (plus Mechanical Check if -1)
495      !!                n = 2    -> exhaustive state
496      !!                n = 3    -> ice/ocean salt fluxes
497      !!
498      !! ** input   :   point coordinates (i,j)
499      !!                n : number of the option
500      !!-------------------------------------------------------------------
501      INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices
502      CHARACTER(len=*), INTENT(in) ::   cd1           !
503      !!
504      INTEGER :: jl
505      !!-------------------------------------------------------------------
506
507      WRITE(numout,*) cd1             ! print title
508
509      !----------------
510      !  Simple state
511      !----------------
512
513      IF ( kn == 1 .OR. kn == -1 ) THEN
514         WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj
515         WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
516         WRITE(numout,*) ' Simple state '
517         WRITE(numout,*) ' masks s,u,v   : ', tms(ki,kj), tmu(ki,kj), tmv(ki,kj)
518         WRITE(numout,*) ' lat - long    : ', gphit(ki,kj), glamt(ki,kj)
519         WRITE(numout,*) ' Time step     : ', numit
520         WRITE(numout,*) ' - Ice drift   '
521         WRITE(numout,*) '   ~~~~~~~~~~~ '
522         WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ki-1,kj)
523         WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ki,kj)
524         WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ki,kj-1)
525         WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ki,kj)
526         WRITE(numout,*) ' strength      : ', strength(ki,kj)
527         WRITE(numout,*)
528         WRITE(numout,*) ' - Cell values '
529         WRITE(numout,*) '   ~~~~~~~~~~~ '
530         WRITE(numout,*) ' cell area     : ', area(ki,kj)
531         WRITE(numout,*) ' at_i          : ', at_i(ki,kj)       
532         WRITE(numout,*) ' vt_i          : ', vt_i(ki,kj)       
533         WRITE(numout,*) ' vt_s          : ', vt_s(ki,kj)       
534         DO jl = 1, jpl
535            WRITE(numout,*) ' - Category (', jl,')'
536            WRITE(numout,*) ' a_i           : ', a_i(ki,kj,jl)
537            WRITE(numout,*) ' ht_i          : ', ht_i(ki,kj,jl)
538            WRITE(numout,*) ' ht_s          : ', ht_s(ki,kj,jl)
539            WRITE(numout,*) ' v_i           : ', v_i(ki,kj,jl)
540            WRITE(numout,*) ' v_s           : ', v_s(ki,kj,jl)
541            WRITE(numout,*) ' e_s           : ', e_s(ki,kj,1,jl)/1.0e9
542            WRITE(numout,*) ' e_i           : ', e_i(ki,kj,1:nlay_i,jl)/1.0e9
543            WRITE(numout,*) ' t_su          : ', t_su(ki,kj,jl)
544            WRITE(numout,*) ' t_snow        : ', t_s(ki,kj,1,jl)
545            WRITE(numout,*) ' t_i           : ', t_i(ki,kj,1:nlay_i,jl)
546            WRITE(numout,*) ' sm_i          : ', sm_i(ki,kj,jl)
547            WRITE(numout,*) ' smv_i         : ', smv_i(ki,kj,jl)
548            WRITE(numout,*)
549            WRITE(numout,*) ' Pathological case : ', patho_case(ki,kj,jl)
550         END DO
551      ENDIF
552      IF( kn == -1 ) THEN
553         WRITE(numout,*) ' Mechanical Check ************** '
554         WRITE(numout,*) ' Check what means ice divergence '
555         WRITE(numout,*) ' Total ice concentration ', at_i (ki,kj)
556         WRITE(numout,*) ' Total lead fraction     ', ato_i(ki,kj)
557         WRITE(numout,*) ' Sum of both             ', ato_i(ki,kj) + at_i(ki,kj)
558         WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ki,kj) + at_i(ki,kj) - 1.00
559      ENDIF
560
561
562      !--------------------
563      !  Exhaustive state
564      !--------------------
565
566      IF ( kn .EQ. 2 ) THEN
567         WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj
568         WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
569         WRITE(numout,*) ' Exhaustive state '
570         WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj)
571         WRITE(numout,*) ' Time step ', numit
572         WRITE(numout,*) 
573         WRITE(numout,*) ' - Cell values '
574         WRITE(numout,*) '   ~~~~~~~~~~~ '
575         WRITE(numout,*) ' cell area     : ', area(ki,kj)
576         WRITE(numout,*) ' at_i          : ', at_i(ki,kj)       
577         WRITE(numout,*) ' vt_i          : ', vt_i(ki,kj)       
578         WRITE(numout,*) ' vt_s          : ', vt_s(ki,kj)       
579         WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ki-1,kj)
580         WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ki,kj)
581         WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ki,kj-1)
582         WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ki,kj)
583         WRITE(numout,*) ' strength      : ', strength(ki,kj)
584         WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ki,kj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ki,kj)
585         WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ki,kj)  , ' old_v_ice     : ', old_v_ice(ki,kj) 
586         WRITE(numout,*)
587
588         DO jl = 1, jpl
589              WRITE(numout,*) ' - Category (',jl,')'
590              WRITE(numout,*) '   ~~~~~~~~         ' 
591              WRITE(numout,*) ' ht_i       : ', ht_i(ki,kj,jl)             , ' ht_s       : ', ht_s(ki,kj,jl)
592              WRITE(numout,*) ' t_i        : ', t_i(ki,kj,1:nlay_i,jl)
593              WRITE(numout,*) ' t_su       : ', t_su(ki,kj,jl)             , ' t_s        : ', t_s(ki,kj,1,jl)
594              WRITE(numout,*) ' sm_i       : ', sm_i(ki,kj,jl)             , ' o_i        : ', o_i(ki,kj,jl)
595              WRITE(numout,*) ' a_i        : ', a_i(ki,kj,jl)              , ' old_a_i    : ', old_a_i(ki,kj,jl)   
596              WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ki,kj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ki,kj,jl) 
597              WRITE(numout,*) ' v_i        : ', v_i(ki,kj,jl)              , ' old_v_i    : ', old_v_i(ki,kj,jl)   
598              WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ki,kj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ki,kj,jl) 
599              WRITE(numout,*) ' v_s        : ', v_s(ki,kj,jl)              , ' old_v_s    : ', old_v_s(ki,kj,jl) 
600              WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ki,kj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ki,kj,jl)
601              WRITE(numout,*) ' e_i1       : ', e_i(ki,kj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ki,kj,1,jl)/1.0e9 
602              WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ki,kj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ki,kj,1,jl)/1.0e9
603              WRITE(numout,*) ' e_i2       : ', e_i(ki,kj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ki,kj,2,jl)/1.0e9 
604              WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ki,kj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ki,kj,2,jl)/1.0e9
605              WRITE(numout,*) ' e_snow     : ', e_s(ki,kj,1,jl)            , ' old_e_snow : ', old_e_s(ki,kj,1,jl) 
606              WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ki,kj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ki,kj,1,jl)
607              WRITE(numout,*) ' smv_i      : ', smv_i(ki,kj,jl)            , ' old_smv_i  : ', old_smv_i(ki,kj,jl)   
608              WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ki,kj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ki,kj,jl) 
609              WRITE(numout,*) ' oa_i       : ', oa_i(ki,kj,jl)             , ' old_oa_i   : ', old_oa_i(ki,kj,jl)
610              WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ki,kj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ki,kj,jl)
611              WRITE(numout,*) ' Path. case : ', patho_case(ki,kj,jl)
612        END DO !jl
613
614        WRITE(numout,*)
615        WRITE(numout,*) ' - Heat / FW fluxes '
616        WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
617!       WRITE(numout,*) ' fsbri      : ', fsbri(ki,kj)
618!       WRITE(numout,*) ' fseqv      : ', fseqv(ki,kj)
619!       WRITE(numout,*) ' fsalt_res  : ', fsalt_res(ki,kj)
620        WRITE(numout,*) ' fmmec      : ', fmmec(ki,kj)
621        WRITE(numout,*) ' fhmec      : ', fhmec(ki,kj)
622        WRITE(numout,*) ' fhbri      : ', fhbri(ki,kj)
623        WRITE(numout,*) ' fheat_rpo  : ', fheat_rpo(ki,kj)
624        WRITE(numout,*) 
625        WRITE(numout,*) ' sst        : ', sst_m(ki,kj) 
626        WRITE(numout,*) ' sss        : ', sss_m(ki,kj) 
627        WRITE(numout,*) 
628        WRITE(numout,*) ' - Stresses '
629        WRITE(numout,*) '   ~~~~~~~~ '
630        WRITE(numout,*) ' utau_ice   : ', utau_ice(ki,kj) 
631        WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ki,kj)
632        WRITE(numout,*) ' utau       : ', utau(ki,kj) 
633        WRITE(numout,*) ' vtau       : ', vtau(ki,kj)
634        WRITE(numout,*) ' oc. vel. u : ', u_oce(ki,kj)
635        WRITE(numout,*) ' oc. vel. v : ', v_oce(ki,kj)
636     ENDIF
637
638     !---------------------
639     ! Salt / heat fluxes
640     !---------------------
641
642     IF ( kn .EQ. 3 ) THEN
643        WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj
644        WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
645        WRITE(numout,*) ' - Salt / Heat Fluxes '
646        WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
647        WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj)
648        WRITE(numout,*) ' Time step ', numit
649        WRITE(numout,*)
650        WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
651        WRITE(numout,*) ' qsr        : ', qsr(ki,kj)
652        WRITE(numout,*) ' qns        : ', qns(ki,kj)
653        WRITE(numout,*)
654        WRITE(numout,*) ' - Salt fluxes at bottom interface ***'
655        WRITE(numout,*) ' emps       : ', emps(ki,kj)
656        WRITE(numout,*) ' emp        : ', emp(ki,kj)
657        WRITE(numout,*) ' fsbri      : ', fsbri(ki,kj)
658        WRITE(numout,*) ' fseqv      : ', fseqv(ki,kj)
659        WRITE(numout,*) ' fsalt_res  : ', fsalt_res(ki,kj)
660        WRITE(numout,*) ' fsalt_rpo  : ', fsalt_rpo(ki,kj)
661        WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
662        WRITE(numout,*) ' fheat_res  : ', fheat_res(ki,kj)
663        WRITE(numout,*)
664        WRITE(numout,*) ' - Momentum fluxes '
665        WRITE(numout,*) ' utau      : ', utau(ki,kj) 
666        WRITE(numout,*) ' vtau      : ', vtau(ki,kj)
667      ENDIF
668      WRITE(numout,*) ' '
669      !
670   END SUBROUTINE lim_prt_state
671
672#else
673   !!----------------------------------------------------------------------
674   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
675   !!----------------------------------------------------------------------
676CONTAINS
677   SUBROUTINE sbc_ice_lim ( kt, kblk )     ! Dummy routine
678      WRITE(*,*) 'sbc_ice_lim: You should not have seen this print! error?', kt, kblk
679   END SUBROUTINE sbc_ice_lim
680#endif
681
682   !!======================================================================
683END MODULE sbcice_lim
Note: See TracBrowser for help on using the repository browser.