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.
icestp.F90 in trunk/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMO/LIM_SRC_3/icestp.F90 @ 833

Last change on this file since 833 was 825, checked in by ctlod, 16 years ago

dev_002_LIM : add the LIM 3.0 component, see ticketr: #71

File size: 37.3 KB
Line 
1MODULE icestp
2   !!======================================================================
3   !!                       ***  MODULE icestp   ***
4   !!   Sea-Ice model : LIM Sea ice model time-stepping
5   !!======================================================================
6#if defined key_lim3
7   !!----------------------------------------------------------------------
8   !!   'key_lim3' :                                   Lim sea-ice model
9   !!----------------------------------------------------------------------
10   !!   ice_stp       : sea-ice model time-stepping
11   !!----------------------------------------------------------------------
12   USE par_ice
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         ! nemo add
20   USE daymod
21   USE phycst          ! Define parameters for the routines
22   USE taumod
23   USE ice
24   USE iceini
25   USE ocesbc
26   USE lbclnk
27   USE limdyn
28   USE limtrp
29   USE limicepoints
30   USE limthd
31   USE limitd_th
32   USE limitd_me
33   USE limflx
34   USE limdia
35   USE limwri
36   USE limrst
37   USE limupdate       ! update of global variables
38   USE limvar
39   USE prtctl          ! Print control
40
41   IMPLICIT NONE
42   PRIVATE
43
44   !! * Routine accessibility
45   PUBLIC ice_stp  ! called by step.F90
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49#  include "vectopt_loop_substitute.h90"
50   !!-----------------------------------------------------
51   !!   LIM 2.0 , UCL-LODYC-IPSL  (2003)
52   !!-----------------------------------------------------
53
54CONTAINS
55
56   SUBROUTINE ice_stp ( kt )
57      !!---------------------------------------------------------------------
58      !!                  ***  ROUTINE ice_stp  ***
59      !!                   
60      !! ** Purpose :   Louvain la Neuve Sea Ice Model time stepping
61      !!
62      !! ** Action  : - call the ice dynamics routine
63      !!              - call the ice advection/diffusion routine
64      !!              - call the ice thermodynamics routine
65      !!              - call the routine that computes mass and
66      !!                heat fluxes at the ice/ocean interface
67      !!              - save the outputs
68      !!              - save the outputs for restart when necessary
69      !!
70      !! History :
71      !!   1.0  !  99-11  (M. Imbard)  Original code
72      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
73      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
74      !!   3.0  !  05-10  (M. Vancoppenolle)   Trend terms, ice salinity, multi-layer model
75      !!----------------------------------------------------------------------
76      !! * Arguments
77      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
78
79      !! * Local declarations
80      INTEGER   ::              & 
81         ji, jj, jk, jl,        &  ! loop indexes
82         index,                 &  ! indexes for ice points
83         numaltests,            &  ! number of alert tests (max 20)
84         alert_id                  ! number of the current alert
85
86      INTEGER, DIMENSION(20) :: &
87         numal                     ! number of alerts positive
88
89      REAL(wp) ::               &
90         ztmelts                   ! ice layer melting point
91      REAL(wp) , DIMENSION(jpi,jpj)    :: &
92         zsss_io, zsss2_io, zsss3_io          ! tempory workspaces
93
94      CHARACTER (len=30), DIMENSION(20) :: &  ! name of alert
95         alname
96     
97      !!----------------------------------------------------------------------
98
99!-------------------------------------------------------------------------------
100! 1) First time step
101!-------------------------------------------------------------------------------
102      !+++++
103      index = 12
104      jiindex = arc_sp_grid(index,1)
105      jjindex = arc_sp_grid(index,2)
106      jiindex = 44
107      jjindex = 140
108
109      WRITE(numout,*) ' The debugging point is : jiindex : ',jiindex,  &
110                                               ' jjindex : ',jjindex
111      !+++++
112
113      IF( kt == nit000 ) THEN
114
115         IF( lk_cpl ) THEN
116            IF(lwp) WRITE(numout,*)
117            IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)'
118            IF(lwp) WRITE(numout,*) '~~~~~~~   coupled case'
119         ELSE
120            IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)' 
121            IF(lwp) WRITE(numout,*) '~~~~~~~   forced case using bulk formulea'
122         ENDIF
123
124         gtaux(:,:) = 0.e0  ! Set wind stress to 0
125         gtauy(:,:) = 0.e0
126
127      ENDIF ! Initial time step
128
129      IF(lwp) WRITE(numout,*) ' ~~~ LIM-@ ~~~      '
130      IF(lwp) WRITE(numout,*) ' Time step : ', numit
131
132!-------------------------------------------------------------------------------
133! 2) SST, SSS, wind stress, ocean velocity and seawater freezing point
134!-------------------------------------------------------------------------------
135
136      !-------------------------------------
137      ! Update SST, SSS and ocean velocity
138      !-------------------------------------
139
140      sst_io(:,:) = sst_io(:,:) + tn(:,:,1) + rt0
141      sss_io(:,:) = sss_io(:,:) + sn(:,:,1)
142
143      ! vectors at F-point
144      DO jj = 2, jpj
145         DO ji = fs_2, jpi 
146            u_io(ji,jj) = u_io(ji,jj) + un(ji,jj,1) 
147            v_io(ji,jj) = v_io(ji,jj) + vn(ji,jj,1) 
148         END DO
149      END DO
150     
151      !------------------------
152      ! Ice model loop starts           
153      !------------------------
154
155      IF ( MOD( kt-1, nfice ) == 0 ) THEN
156         
157         !-------------------
158         ! Average SST, SSS
159         !-------------------
160
161         sst_io(:,:) = sst_io(:,:) / FLOAT( nfice ) * tmask(:,:,1)
162         sss_io(:,:) = sss_io(:,:) / FLOAT( nfice )
163
164         !----------------------------------------
165         ! Atmospheric stress and ocean velocity
166         !----------------------------------------
167
168         DO jj = 2, jpj
169            DO ji = fs_2, jpi
170
171               gtaux(ji,jj) = cai / cao * tauxw(ji,jj)
172               gtauy(ji,jj) = cai / cao * tauyw(ji,jj)
173
174               u_io (ji,jj) = u_io(ji,jj) / FLOAT( nfice ) * tmu(ji,jj)
175               v_io (ji,jj) = v_io(ji,jj) / FLOAT( nfice ) * tmv(ji,jj)
176
177            END DO
178         END DO
179
180         ! Boundary conditions
181         CALL lbc_lnk( gtaux(:,:), 'U', -1. )   ! I-point (i.e. ice U-V point)
182         CALL lbc_lnk( gtauy(:,:), 'V', -1. )   ! I-point (i.e. ice U-V point)
183         CALL lbc_lnk( u_io (:,:), 'U', -1. )   ! I-point (i.e. ice U-V point)
184         CALL lbc_lnk( v_io (:,:), 'V', -1. )   ! I-point (i.e. ice U-V point)
185
186         !--------------------------
187         ! Seawater freezing point
188         !--------------------------
189
190         zsss_io (:,:) = SQRT( sss_io(:,:) ) 
191         zsss2_io(:,:) = sss_io(:,:) * sss_io(:,:)
192         zsss3_io(:,:) = zsss_io(:,:) * zsss_io(:,:) * zsss_io(:,:)
193
194         DO jj = 1, jpj
195            DO ji = 1, jpi
196               t_bo(ji,jj)  = ABS ( rt0 - 0.0575      * sss_io(ji,jj)      &
197                  &                    + 1.710523e-03 * zsss3_io(ji,jj)    &
198                  &                    - 2.154996e-04 * zsss2_io(ji,jj) )  &
199                            * tms(ji,jj)
200            END DO
201         END DO
202
203         !+++++
204         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
205            CALL prt_ctl_info('Ice Forcings ')
206            CALL prt_ctl(tab2d_1=qsr_oce ,clinfo1=' qsr_oce  : ')
207            CALL prt_ctl(tab2d_1=qnsr_oce,clinfo1=' qnsr_oce : ')
208            CALL prt_ctl(tab2d_1=evap    ,clinfo1=' evap     : ')
209            CALL prt_ctl(tab2d_1=tprecip ,clinfo1=' precip   : ')
210            CALL prt_ctl(tab2d_1=gtaux   ,clinfo1=' u-stress : ', tab2d_2=gtauy   , clinfo2=' v-stress  : ')
211            CALL prt_ctl(tab2d_1=sst_io  ,clinfo1=' sst      : ', tab2d_2=sss_io  , clinfo2=' sss       : ')
212            CALL prt_ctl(tab2d_1=u_io    ,clinfo1=' u_io     : ', tab2d_2=v_io    , clinfo2=' v_io      : ')
213            CALL prt_ctl(tab2d_1=frld    ,clinfo1=' frld   1 : ')
214
215            DO jl = 1, jpl
216               CALL prt_ctl_info('* - category number ', ivar1=jl)
217               CALL prt_ctl(tab3d_1=t_su    , clinfo1=' t_su      : ', kdim=jl)
218               CALL prt_ctl(tab3d_1=qsr_ice , clinfo1=' qsr_ice   : ', kdim=jl)
219               CALL prt_ctl(tab3d_1=qnsr_ice, clinfo1=' qnsr_ice  : ', kdim=jl)
220               CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' dqns_ice  : ', kdim=jl)
221               CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' qla_ice   : ', kdim=jl)
222               CALL prt_ctl(tab3d_1=dqla_ice, clinfo1=' dqla_ice  : ', kdim=jl)
223            END DO
224         ENDIF
225         !+++++
226
227!----------------------------------------------------------------------------
228! 3) Store old values of ice model global variables
229!----------------------------------------------------------------------------
230
231         old_a_i(:,:,:)   = a_i(:,:,:)     ! ice area
232         old_e_i(:,:,:,:) = e_i(:,:,:,:)   ! ice thermal energy
233         old_v_i(:,:,:)   = v_i(:,:,:)     ! ice volume
234         old_v_s(:,:,:)   = v_s(:,:,:)     ! snow volume
235         old_e_s(:,:,:,:) = e_s(:,:,:,:)   ! snow thermal energy
236         old_smv_i(:,:,:) = smv_i(:,:,:)   ! salt content
237         old_oa_i(:,:,:)  = oa_i(:,:,:)    ! areal age content
238
239         d_a_i_thd(:,:,:)   = 0.0 ; d_a_i_trp(:,:,:)   = 0.0
240         d_v_i_thd(:,:,:)   = 0.0 ; d_v_i_trp(:,:,:)   = 0.0
241         d_e_i_thd(:,:,:,:) = 0.0 ; d_e_i_trp(:,:,:,:) = 0.0
242         d_v_s_thd(:,:,:)   = 0.0 ; d_v_s_trp(:,:,:)   = 0.0
243         d_e_s_thd(:,:,:,:) = 0.0 ; d_e_s_trp(:,:,:,:) = 0.0
244         d_smv_i_thd(:,:,:) = 0.0 ; d_smv_i_trp(:,:,:) = 0.0
245         d_oa_i_thd(:,:,:)  = 0.0 ; d_oa_i_trp(:,:,:)  = 0.0
246
247         fsalt(:,:)     = 0.0     ; fseqv(:,:)     = 0.0
248         fsbri(:,:)     = 0.0     ; fsalt_res(:,:) = 0.0
249         fsalt_rpo(:,:) = 0.0
250         fhmec(:,:)     = 0.0     ; fhbri(:,:) = 0.0; 
251         fmmec(:,:)     = 0.0     ; fheat_res(:,:) = 0.0
252         fheat_rpo(:,:) = 0.0     ; focea2D(:,:) = 0.0
253         fsup2D(:,:)    = 0.0
254
255         diag_sni_gr(:,:) = 0.0   ; diag_lat_gr(:,:) = 0.0
256         diag_bot_gr(:,:) = 0.0   ; diag_dyn_gr(:,:) = 0.0
257         diag_bot_me(:,:) = 0.0   ; diag_sur_me(:,:) = 0.0
258
259         ! dynamical invariants
260         delta_i(:,:) = 0.0
261         divu_i(:,:)  = 0.0
262         shear_i(:,:) = 0.0
263
264         !----------------------
265         ! Ice model time step
266         !----------------------
267         numit = numit + nfice 
268         
269!----------------------------------------------------------------------------
270! 4) Ice routines, compute the trend terms of global variables
271!----------------------------------------------------------------------------
272
273         !+++++
274         WRITE(numout,*) ' - Beginning the time step - '
275         CALL lim_inst_state(jiindex,jjindex,1)
276         WRITE(numout,*) ' '
277         !+++++
278
279        !---------------------------
280        ! 4.1) Dynamical processes         
281        !---------------------------
282
283        !                            !--------------------!
284        CALL lim_dyn                 !  Ice dynamics      !  ( rheology/dynamics )
285        !                            !--------------------!
286
287        !                            !--------------------!
288        CALL lim_trp                 !  Ice transport     !
289        !                            !--------------------!
290   
291        CALL lim_var_agg(1)  ! aggregate categories, requested
292        CALL lim_var_glo2eqv ! equivalent variables, requested for rafting
293
294        !+++++
295        WRITE(numout,*) ' - After ice dynamics and transport '
296        CALL lim_inst_state(jiindex,jjindex,1)
297        WRITE(numout,*) ' '
298        WRITE(numout,*)
299        WRITE(numout,*) ' Mechanical Check ************** '
300        WRITE(numout,*) ' Check what means ice divergence '
301        WRITE(numout,*) ' Total ice concentration ', at_i (jiindex,jjindex)
302        WRITE(numout,*) ' Total lead fraction     ', ato_i(jiindex,jjindex)
303        WRITE(numout,*) ' Sum of both             ', ato_i(jiindex,jjindex) + at_i(jiindex,jjindex)
304        WRITE(numout,*) ' Sum of both minus 1     ', ato_i(jiindex,jjindex) + at_i(jiindex,jjindex) - 1.00
305        !+++++
306
307        !                            !------------------------------------------------!
308        CALL lim_itd_me              !  Mechanical redistribution ! (ridging/rafting)
309        !                            !------------------------------------------------!
310
311        !-------------------------
312        ! 4.2 Ice thermodynamics
313        !-------------------------
314
315        CALL lim_var_glo2eqv ! equivalent variables
316        CALL lim_var_agg(1)  ! aggregate ice categories
317        CALL lim_var_bv      ! bulk brine volume (diag)
318
319        !                            !--------------------------------!
320        CALL lim_thd                 !  Heat diffusion, growth, melt  !
321        !                            !--------------------------------!
322     
323        !                            !---------------------!
324        !                            !  Ice natural aging  !
325        !                            !---------------------!
326
327        oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice / 86400.00
328
329        !+++++
330        CALL lim_var_glo2eqv ! except the info message that follows, i dont
331                             ! think this one is necessary
332        WRITE(numout,*) ' - After ice thermodynamics '
333        CALL lim_inst_state(jiindex,jjindex,1)
334        !+++++
335
336
337        !                            !-------------------------------------------!
338        CALL lim_itd_th              !  Remap ice categories, lateral accretion  !
339        !                            !-------------------------------------------!
340
341!----------------------------------------------------------------------------
342! 5) Global variables update
343!----------------------------------------------------------------------------
344
345        CALL lim_var_agg(1)   ! requested by limupdate
346        !                            !-------------------------!
347        CALL lim_update              ! Global variables update !
348        !                            !-------------------------!
349
350        CALL lim_var_glo2eqv ! equivalent variables (outputs)
351        CALL lim_var_agg(2)   ! aggregate ice thickness categories
352
353        !+++++
354        IF(ln_nicep) THEN
355           WRITE(numout,*) ' - Final ice state after lim_update '
356           CALL lim_inst_state(jiindex,jjindex,2)
357           WRITE(numout,*) ' '
358        ENDIF
359        !+++++
360
361!----------------------------------------------------------------------------
362! 6) Fluxes of mass and heat to the ocean
363!----------------------------------------------------------------------------
364
365        !                        !------------------------------!
366        CALL lim_flx             ! Ice/Ocean Mass & Heat fluxes !
367        !                        !------------------------------!
368
369        !+++++
370        WRITE(numout,*) ' - Final ice state after lim_flx    '
371        CALL lim_inst_state(jiindex,jjindex,3)
372        WRITE(numout,*) ' '
373        !+++++
374
375!----------------------------------------------------------------------------
376! 7) Diagnostics and outputs
377!----------------------------------------------------------------------------
378
379         IF( MOD( numit, ninfo ) == 0 .OR. ntmoy == 1 )  THEN       
380         !                        !-------------------------------!
381         CALL lim_dia             ! Ice Diagnostics in evolu file !
382         ENDIF                    !-------------------------------!
383
384         !                        !----------------------------!
385         CALL lim_wri( 1 )        ! Ice outputs in icemod file !
386         !                        !----------------------------!
387
388         IF( MOD( numit, nstock ) == 0 .OR. numit == nlast ) THEN
389                 WRITE(numout,*) ' Ice restart is called '
390         !                           !------------------!
391         CALL lim_rst_write( numit ) ! Ice restart file !
392         !                           !------------------!
393         ENDIF
394
395         CALL lim_var_glo2eqv
396
397         ! Re-initialization of forcings
398         qsr_oce (:,:) = 0.e0
399         qsr_ice (:,:,:) = 0.e0
400         qnsr_oce(:,:) = 0.e0
401         qnsr_ice(:,:,:) = 0.e0 
402         dqns_ice(:,:,:) = 0.e0 
403         tprecip (:,:) = 0.e0 
404         sprecip (:,:) = 0.e0
405#if defined key_coupled 
406         rrunoff (:,:) = 0.e0
407         calving (:,:) = 0.e0
408#else
409         qla_ice (:,:,:) = 0.e0
410         dqla_ice(:,:,:) = 0.e0
411         fr1_i0  (:,:) = 0.e0
412         fr2_i0  (:,:) = 0.e0
413         evap    (:,:) = 0.e0
414#endif
415
416!----------------------------------------------------------------------------
417! 7) Alerts in case of model crash
418!----------------------------------------------------------------------------
419
420     numaltests = 10
421     numal(:) = 0
422
423     ! Alert if incompatible volume and concentration
424     alert_id = 2 ! reference number of this alert
425     alname(alert_id) = ' Incompat vol and con         ' ! name of the alert
426
427     DO jl = 1, jpl
428        DO jj = 1, jpj
429           DO ji = 1, jpi
430                 IF ((v_i(ji,jj,jl).NE.0.0).AND.(a_i(ji,jj,jl).EQ.0.0)) THEN
431                    WRITE(numout,*) ' ALERTE 2 '
432                    WRITE(numout,*) ' Incompatible volume and concentration '
433                    WRITE(numout,*) ' at_i     ', at_i(ji,jj)
434                    WRITE(numout,*) ' Point - category', ji, jj, jl
435                    WRITE(numout,*)
436                    WRITE(numout,*) ' a_i *** a_i_old ', a_i(ji,jj,jl), old_a_i(ji,jj,jl)
437                    WRITE(numout,*) ' v_i *** v_i_old ', v_i(ji,jj,jl), old_v_i(ji,jj,jl)
438                    WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
439                    WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
440                    numal(alert_id) = numal(alert_id) + 1
441                 ENDIF
442           END DO
443        END DO
444     END DO
445
446     ! Alerte if very thick ice
447     alert_id = 3 ! reference number of this alert
448     alname(alert_id) = ' Very thick ice               ' ! name of the alert
449     jl = jpl 
450     DO jj = 1, jpj
451        DO ji = 1, jpi
452           IF (ht_i(ji,jj,jl) .GT. 50.0 ) THEN
453              WRITE(numout,*) ' ALERTE 3 '
454              WRITE(numout,*) ' Very thick ice '
455              CALL lim_inst_state(ji,jj,2)
456              WRITE(numout,*) ' '
457              numal(alert_id) = numal(alert_id) + 1
458           ENDIF
459        END DO
460     END DO
461
462     ! Alert if very fast ice
463     alert_id = 4 ! reference number of this alert
464     alname(alert_id) = ' Very fast ice               ' ! name of the alert
465     DO jj = 1, jpj
466        DO ji = 1, jpi
467           IF ( ( ( ABS( u_ice(ji,jj) ) .GT. 0.50) .OR. &
468                  ( ABS( v_ice(ji,jj) ) .GT. 0.50) ) .AND. &
469                ( at_i(ji,jj) .GT. 0.0 ) ) THEN
470              WRITE(numout,*) ' ALERTE 4 '
471              WRITE(numout,*) ' Very fast ice '
472              CALL lim_inst_state(ji,jj,1)
473              WRITE(numout,*) ' ice strength         : ', strength(ji,jj)
474              WRITE(numout,*) ' oceanic stress ftaux : ', ftaux(ji,jj) 
475              WRITE(numout,*) ' oceanic stress ftauy : ', ftauy(ji,jj)
476              WRITE(numout,*) ' atmosph stress gtaux : ', gtaux(ji,jj) 
477              WRITE(numout,*) ' atmosph stress gtauy : ', gtauy(ji,jj)
478              WRITE(numout,*) ' oceanic speed u      : ', u_io(ji,jj)
479              WRITE(numout,*) ' oceanic speed v      : ', v_io(ji,jj)
480              WRITE(numout,*) ' sst                  : ', sst_io(ji,jj)
481              WRITE(numout,*) ' sss                  : ', sss_io(ji,jj)
482              WRITE(numout,*) 
483              numal(alert_id) = numal(alert_id) + 1
484           ENDIF
485        END DO
486     END DO
487
488     ! Alert if there is ice on continents
489     alert_id = 6 ! reference number of this alert
490     alname(alert_id) = ' Ice on continents           ' ! name of the alert
491     DO jj = 1, jpj
492        DO ji = 1, jpi
493           IF ( ( tms(ji,jj) .LE. 0.0 ) .AND. ( at_i(ji,jj) .GT. 0.0 ) ) THEN
494              WRITE(numout,*) ' ALERTE 6 '
495              WRITE(numout,*) ' Ice on continents '
496              CALL lim_inst_state(ji,jj,1)
497              WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), &
498                                                          tmu(ji,jj), &
499                                                          tmv(ji,jj) 
500              WRITE(numout,*) ' sst                  : ', sst_io(ji,jj)
501              WRITE(numout,*) ' sss                  : ', sss_io(ji,jj)
502              WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj)
503              WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj)
504              WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1)
505              WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj)
506              WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj)
507
508              numal(alert_id) = numal(alert_id) + 1
509
510           ENDIF
511        END DO
512     END DO
513
514!
515!    ! Alert if very fresh ice
516     alert_id = 7 ! reference number of this alert
517     alname(alert_id) = ' Very fresh ice               ' ! name of the alert
518     DO jl = 1, jpl
519     DO jj = 1, jpj
520        DO ji = 1, jpi
521           IF ( ( ( ABS( sm_i(ji,jj,jl) ) .LT. 0.50) .OR. &
522                  ( ABS( sm_i(ji,jj,jl) ) .LT. 0.50) ) .AND. &
523                ( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN
524!             WRITE(numout,*) ' ALERTE 7 '
525!             WRITE(numout,*) ' Very fresh ice '
526!             CALL lim_inst_state(ji,jj,1)
527!             WRITE(numout,*) ' sst                  : ', sst_io(ji,jj)
528!             WRITE(numout,*) ' sss                  : ', sss_io(ji,jj)
529!             WRITE(numout,*) ' s_i_newice           : ', s_i_newice(ji,jj,1:jpl)
530!             WRITE(numout,*)
531              numal(alert_id) = numal(alert_id) + 1
532           ENDIF
533        END DO
534     END DO
535     END DO
536!
537
538!    ! Alert if too old ice
539     alert_id = 9 ! reference number of this alert
540     alname(alert_id) = ' Very old   ice               ' ! name of the alert
541     DO jl = 1, jpl
542     DO jj = 1, jpj
543        DO ji = 1, jpi
544           IF ( ( ( ABS( o_i(ji,jj,jl) ) .GT. rdt_ice ) .OR. &
545                  ( ABS( o_i(ji,jj,jl) ) .LT. 0.00) ) .AND. &
546                ( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN
547              WRITE(numout,*) ' ALERTE 9 '
548              WRITE(numout,*) ' Wrong ice age '
549              CALL lim_inst_state(ji,jj,1)
550              WRITE(numout,*) 
551              numal(alert_id) = numal(alert_id) + 1
552           ENDIF
553        END DO
554     END DO
555
556     END DO
557
558     ! Alert on salt flux
559     alert_id = 5 ! reference number of this alert
560     alname(alert_id) = ' High salt flux               ' ! name of the alert
561     DO jj = 1, jpj
562        DO ji = 1, jpi
563           IF (ABS(fsalt(ji,jj)).gt.1.0e-2) THEN
564           WRITE(numout,*) ' ALERTE 5 '
565           WRITE(numout,*) ' High salt flux '
566           CALL lim_inst_state(ji,jj,3)
567           WRITE(numout,*) ' '
568           DO jl = 1, jpl
569              WRITE(numout,*) ' Category no: ', jl
570              WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)     , &
571                              ' old_a_i    : ', old_a_i(ji,jj,jl)   
572              WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , &
573                              ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl) 
574              WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)    , &
575                              ' old_v_i    : ', old_v_i(ji,jj,jl)   
576              WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , &
577                              ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl) 
578              WRITE(numout,*) ' '
579              END DO
580              numal(alert_id) = numal(alert_id) + 1
581           ENDIF
582        END DO
583     END DO
584
585     ! Alert if fnsolar very big
586     alert_id = 8 ! reference number of this alert
587     alname(alert_id) = ' fnsolar very big             ' ! name of the alert
588     DO jj = 1, jpj
589        DO ji = 1, jpi
590           IF ( ( ABS( fnsolar(ji,jj) ) .GT. 1500.0 ) .AND. &
591                ( at_i(ji,jj) .GT. 0.0 ) )  THEN
592
593              WRITE(numout,*) ' ALERTE 8 '
594              WRITE(numout,*) ' ji, jj    : ', ji, jj
595              WRITE(numout,*) ' fnsolar : ', fnsolar(ji,jj)
596              WRITE(numout,*) ' Very high non-solar heat flux '
597              WRITE(numout,*) ' sst       : ', sst_io(ji,jj)
598              WRITE(numout,*) ' sss       : ', sss_io(ji,jj)
599              WRITE(numout,*) ' qcmif     : ', qcmif(jiindex,jjindex)
600              WRITE(numout,*) ' qldif     : ', qldif(jiindex,jjindex)
601              WRITE(numout,*) ' qcmif     : ', qcmif(jiindex,jjindex) / rdt_ice
602              WRITE(numout,*) ' qldif     : ', qldif(jiindex,jjindex) / rdt_ice
603              WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindex,jjindex)
604              WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindex,jjindex)
605              WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindex,jjindex) / rdt_ice
606              WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindex,jjindex) / rdt_ice
607              WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindex,jjindex) 
608              WRITE(numout,*) ' fhmec     : ', fhmec(jiindex,jjindex) 
609              WRITE(numout,*) ' fheat_rpo : ', fheat_rpo(jiindex,jjindex) 
610              WRITE(numout,*) ' fheat_res : ', fheat_res(jiindex,jjindex) 
611              WRITE(numout,*) ' fhbri     : ', fhbri(jiindex,jjindex) 
612
613              CALL lim_inst_state(ji,jj,2)
614              numal(alert_id) = numal(alert_id) + 1
615
616           ENDIF
617        END DO
618     END DO
619     !+++++
620 
621     ! Alert if very warm ice
622     alert_id = 10 ! reference number of this alert
623     alname(alert_id) = ' Very warm ice                ' ! name of the alert
624     numal(alert_id) = 0
625     DO jl = 1, jpl
626     DO jk = 1, nlay_i
627     DO jj = 1, jpj
628        DO ji = 1, jpi
629           ztmelts    =  -tmut*s_i(ji,jj,jk,jl) + rtt
630           IF ( ( t_i(ji,jj,jk,jl) .GE. ztmelts) .AND. &
631                ( v_i(ji,jj,jl)   .GT. 1.0e-6)   .AND. &
632                ( a_i(ji,jj,jl) .GT. 0.0     )  ) THEN
633              WRITE(numout,*) ' ALERTE 10 '
634              WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl
635              WRITE(numout,*) ' Very warm ice '
636              WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl)
637              WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
638              WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl)
639              WRITE(numout,*) ' ztmelts : ', ztmelts
640              numal(alert_id) = numal(alert_id) + 1
641           ENDIF
642        END DO
643     END DO
644     END DO
645     END DO
646
647     alert_id = 1 ! reference number of this alert
648     alname(alert_id) = ' Il n''y a pas d''alerte 1       ' ! name of the alert
649     WRITE(numout,*)
650     WRITE(numout,*) ' All alerts at the end of ice model '
651     DO alert_id = 1, numaltests
652        WRITE(numout,*) alert_id, alname(alert_id)//' : ', numal(alert_id), ' times ! '
653     END DO
654!
655
656      ENDIF ! kt EQ ice time step
657
658   END SUBROUTINE ice_stp
659   
660!------------------------------------------------------------------------------
661   SUBROUTINE lim_inst_state(i,j,n)
662      !!-----------------------------------------------------------------------
663      !!                   ***  ROUTINE lim_inst_state ***
664      !!                 
665      !! ** Purpose :   Writes global ice state on the (i,j) point
666      !!                in ocean.ouput
667      !!                3 possibilities exist
668      !!                n = 1 -> simple ice state
669      !!                n = 2 -> exhaustive state
670      !!                n = 3 -> ice/ocean salt fluxes
671      !!
672      !! ** input   :   point coordinates (i,j)
673      !!                n : number of the option
674      !!
675      !! history :
676      !!  9.0  ! 06-12 (M. Vancoppenolle) Original code
677      !!-------------------------------------------------------------------
678      INTEGER, INTENT( in ) ::   i,j,n     ! ocean time-step index
679
680      INTEGER :: jl
681
682      !!-------------------------------------------------------------------
683
684      !----------------
685      !  Simple state
686      !----------------
687
688      IF ( n .EQ. 1 ) THEN
689         WRITE(numout,*) ' lim_inst_state - Point : ',i,j
690         WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
691         WRITE(numout,*) ' Simple state '
692         WRITE(numout,*) ' masks s,u,v   : ', tms(i,j), tmu(i,j), tmv(i,j)
693         WRITE(numout,*) ' lat - long    : ', gphit(i,j), glamt(i,j)
694         WRITE(numout,*) ' Time step     : ', numit
695         WRITE(numout,*) ' - Ice drift   '
696         WRITE(numout,*) '   ~~~~~~~~~~~ '
697         WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(i-1,j)
698         WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(i,j)
699         WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(i,j-1)
700         WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(i,j)
701         WRITE(numout,*) ' strength      : ', strength(i,j)
702         WRITE(numout,*)
703         WRITE(numout,*) ' - Cell values '
704         WRITE(numout,*) '   ~~~~~~~~~~~ '
705         WRITE(numout,*) ' cell area     : ', area(i,j)
706         WRITE(numout,*) ' at_i          : ', at_i(i,j)       
707         WRITE(numout,*) ' vt_i          : ', vt_i(i,j)       
708         WRITE(numout,*) ' vt_s          : ', vt_s(i,j)       
709         DO jl = 1, jpl
710            WRITE(numout,*) ' - Category (', jl,')'
711            WRITE(numout,*) ' a_i           : ', a_i(i,j,jl)
712            WRITE(numout,*) ' ht_i          : ', ht_i(i,j,jl)
713            WRITE(numout,*) ' ht_s          : ', ht_s(i,j,jl)
714            WRITE(numout,*) ' v_i           : ', v_i(i,j,jl)
715            WRITE(numout,*) ' v_s           : ', v_s(i,j,jl)
716            WRITE(numout,*) ' e_s           : ', e_s(i,j,1,jl)/1.0e9
717            WRITE(numout,*) ' e_i           : ', e_i(i,j,1:nlay_i,jl)/1.0e9
718            WRITE(numout,*) ' t_su          : ', t_su(i,j,jl)
719            WRITE(numout,*) ' t_snow        : ', t_s(i,j,1,jl)
720            WRITE(numout,*) ' t_i           : ', t_i(i,j,1:nlay_i,jl)
721            WRITE(numout,*) ' sm_i          : ', sm_i(i,j,jl)
722            WRITE(numout,*) ' smv_i         : ', smv_i(i,j,jl)
723            WRITE(numout,*)
724            WRITE(numout,*) ' Pathological case : ', patho_case(i,j,jl)
725         END DO
726      ENDIF
727
728      !--------------------
729      !  Exhaustive state
730      !--------------------
731
732      IF ( n .EQ. 2 ) THEN
733         WRITE(numout,*) ' lim_inst_state - Point : ',i,j
734         WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
735         WRITE(numout,*) ' Exhaustive state '
736         WRITE(numout,*) ' lat - long ', gphit(i,j), glamt(i,j)
737         WRITE(numout,*) ' Time step ', numit
738         WRITE(numout,*) 
739         WRITE(numout,*) ' - Cell values '
740         WRITE(numout,*) '   ~~~~~~~~~~~ '
741         WRITE(numout,*) ' cell area     : ', area(i,j)
742         WRITE(numout,*) ' at_i          : ', at_i(i,j)       
743         WRITE(numout,*) ' vt_i          : ', vt_i(i,j)       
744         WRITE(numout,*) ' vt_s          : ', vt_s(i,j)       
745         WRITE(numout,*) ' sdvt          : ', sdvt(i,j)       
746         WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(i-1,j)
747         WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(i,j)
748         WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(i,j-1)
749         WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(i,j)
750         WRITE(numout,*) ' strength      : ', strength(i,j)
751         WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(i,j), &
752                         ' d_v_ice_dyn   : ', d_v_ice_dyn(i,j)
753         WRITE(numout,*) ' old_u_ice     : ', old_u_ice(i,j)  , &
754                         ' old_v_ice     : ', old_v_ice(i,j) 
755         WRITE(numout,*)
756
757         DO jl = 1, jpl
758              WRITE(numout,*) ' - Category (',jl,')'
759              WRITE(numout,*) '   ~~~~~~~~         ' 
760              WRITE(numout,*) ' ht_i       : ', ht_i(i,j,jl)             , &
761                              ' ht_s       : ', ht_s(i,j,jl)
762              WRITE(numout,*) ' t_i        : ', t_i(i,j,1:nlay_i,jl)
763              WRITE(numout,*) ' t_su       : ', t_su(i,j,jl)             , &
764                              ' t_s        : ', t_s(i,j,1,jl)
765              WRITE(numout,*) ' sm_i       : ', sm_i(i,j,jl)             , &   
766                              ' o_i        : ', o_i(i,j,jl)
767              WRITE(numout,*) ' a_i        : ', a_i(i,j,jl)              , &
768                              ' old_a_i    : ', old_a_i(i,j,jl)   
769              WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(i,j,jl)        , &
770                              ' d_a_i_thd  : ', d_a_i_thd(i,j,jl) 
771              WRITE(numout,*) ' v_i        : ', v_i(i,j,jl)              , &
772                              ' old_v_i    : ', old_v_i(i,j,jl)   
773              WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(i,j,jl)        , &
774                              ' d_v_i_thd  : ', d_v_i_thd(i,j,jl) 
775              WRITE(numout,*) ' v_s        : ', v_s(i,j,jl)              , &
776                              ' old_v_s    : ', old_v_s(i,j,jl) 
777              WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(i,j,jl)        , &
778                              ' d_v_s_thd  : ', d_v_s_thd(i,j,jl)
779              WRITE(numout,*) ' e_i1       : ', e_i(i,j,1,jl)/1.0e9      , &
780                              ' old_ei1    : ', old_e_i(i,j,1,jl)/1.0e9 
781              WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(i,j,1,jl)/1.0e9, &
782                              ' de_i1_thd  : ', d_e_i_thd(i,j,1,jl)/1.0e9
783              WRITE(numout,*) ' e_i2       : ', e_i(i,j,2,jl)/1.0e9      , &
784                              ' old_ei2    : ', old_e_i(i,j,2,jl)/1.0e9 
785              WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(i,j,2,jl)/1.0e9, &
786                              ' de_i2_thd  : ', d_e_i_thd(i,j,2,jl)/1.0e9
787              WRITE(numout,*) ' e_snow     : ', e_s(i,j,1,jl)            , &
788                              ' old_e_snow : ', old_e_s(i,j,1,jl) 
789              WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(i,j,1,jl)      , &
790                              ' d_e_s_thd  : ', d_e_s_thd(i,j,1,jl)
791              WRITE(numout,*) ' smv_i      : ', smv_i(i,j,jl)            , &
792                              ' old_smv_i  : ', old_smv_i(i,j,jl)   
793              WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(i,j,jl)      , &
794                              ' d_smv_i_thd: ', d_smv_i_thd(i,j,jl) 
795              WRITE(numout,*) ' oa_i       : ', oa_i(i,j,jl)             , &
796                              ' old_oa_i   : ', old_oa_i(i,j,jl)
797              WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(i,j,jl)       , &
798                              ' d_oa_i_thd : ', d_oa_i_thd(i,j,jl)
799              WRITE(numout,*) ' Path. case : ', patho_case(i,j,jl)
800        END DO !jl
801
802        WRITE(numout,*)
803        WRITE(numout,*) ' - Heat / FW fluxes '
804        WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
805!       WRITE(numout,*) ' fsalt      : ', fsalt(i,j)
806!       WRITE(numout,*) ' fmass      : ', fmass(i,j)
807!       WRITE(numout,*) ' fsbri      : ', fsbri(i,j)
808!       WRITE(numout,*) ' fseqv      : ', fseqv(i,j)
809!       WRITE(numout,*) ' fsalt_res  : ', fsalt_res(i,j)
810        WRITE(numout,*) ' fmmec      : ', fmmec(i,j)
811        WRITE(numout,*) ' fhmec      : ', fhmec(i,j)
812        WRITE(numout,*) ' fhbri      : ', fhbri(i,j)
813        WRITE(numout,*) ' fheat_rpo  : ', fheat_rpo(i,j)
814        WRITE(numout,*) 
815        WRITE(numout,*) ' sst        : ', sst_io(i,j) 
816        WRITE(numout,*) ' sss        : ', sss_io(i,j) 
817        WRITE(numout,*) 
818        WRITE(numout,*) ' - Stresses '
819        WRITE(numout,*) '   ~~~~~~~~ '
820        WRITE(numout,*) ' tauxw      : ', tauxw(i,j) 
821        WRITE(numout,*) ' tauyw      : ', tauyw(i,j)
822        WRITE(numout,*) ' taux       : ', taux(i,j) 
823        WRITE(numout,*) ' tauy       : ', tauy(i,j)
824        WRITE(numout,*) ' ftaux      : ', ftaux(i,j) 
825        WRITE(numout,*) ' ftauy      : ', ftauy(i,j)
826        WRITE(numout,*) ' gtaux      : ', gtaux(i,j) 
827        WRITE(numout,*) ' gtauy      : ', gtauy(i,j)
828        WRITE(numout,*) ' oc. vel. u : ', u_io(i,j)
829        WRITE(numout,*) ' oc. vel. v : ', v_io(i,j)
830     ENDIF
831
832     !---------------------
833     ! Salt / heat fluxes
834     !---------------------
835
836     IF ( n .EQ. 3 ) THEN
837        WRITE(numout,*) ' lim_inst_state - Point : ',i,j
838        WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
839        WRITE(numout,*) ' - Salt / Heat Fluxes '
840        WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
841        WRITE(numout,*) ' lat - long ', gphit(i,j), glamt(i,j)
842        WRITE(numout,*) ' Time step ', numit
843        WRITE(numout,*)
844        WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
845        WRITE(numout,*) ' fsolar     : ', fsolar(i,j)
846        WRITE(numout,*) ' fnsolar    : ', fnsolar(i,j)
847        WRITE(numout,*)
848        WRITE(numout,*) ' - Salt fluxes at bottom interface ***'
849        WRITE(numout,*) ' fsalt      : ', fsalt(i,j)
850        WRITE(numout,*) ' fmass      : ', fmass(i,j)
851        WRITE(numout,*) ' fsbri      : ', fsbri(i,j)
852        WRITE(numout,*) ' fseqv      : ', fseqv(i,j)
853        WRITE(numout,*) ' fsalt_res  : ', fsalt_res(i,j)
854        WRITE(numout,*) ' fsalt_rpo  : ', fsalt_rpo(i,j)
855        WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
856        WRITE(numout,*) ' fheat_res  : ', fheat_res(i,j)
857        WRITE(numout,*)
858        WRITE(numout,*) ' - Momentum fluxes '
859        WRITE(numout,*) ' ftaux      : ', ftaux(i,j) 
860        WRITE(numout,*) ' ftauy      : ', ftauy(i,j)
861        WRITE(numout,*) ' gtaux      : ', gtaux(i,j) 
862        WRITE(numout,*) ' gtauy      : ', gtauy(i,j)
863      ENDIF
864
865   END SUBROUTINE
866
867#else
868   !!----------------------------------------------------------------------
869   !!   Default option                                 NO LIM sea-ice model
870   !!----------------------------------------------------------------------
871   USE in_out_manager
872
873CONTAINS
874
875   SUBROUTINE ice_stp ( kt )            ! Empty routine
876      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
877
878      IF( kt == nit000 ) THEN
879         IF(lwp) WRITE(numout,*)
880         IF(lwp) WRITE(numout,*) 'No Sea Ice Model'
881         IF(lwp) WRITE(numout,*) '~~~~~~~'
882      ENDIF
883
884   END SUBROUTINE ice_stp
885   
886   SUBROUTINE lim_inst_state(i,j,n)        ! Empty routine
887      INTEGER, INTENT( in ) :: i, j, n
888   END SUBROUTINE
889
890#endif
891
892   !!======================================================================
893END MODULE icestp
Note: See TracBrowser for help on using the repository browser.