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 @ 835

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

Use IOM interface to build sea-ice restart files, see ticket:#74

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