source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90 @ 8522

Last change on this file since 8522 was 8522, checked in by clem, 3 years ago

changes in style - part6 - ice diffusion (hope the scheme still converges)

File size: 16.2 KB
Line 
1MODULE icedia
2   !!======================================================================
3   !!                       ***  MODULE icedia  ***
4   !!  Sea-Ice model :   global budgets
5   !!======================================================================
6   !! History :  3.4  ! 2012-10  (C. Rousset)  original code
7   !!            4.0  ! 2017-08  (C. Rousset)  fits nemo4.0 standards
8   !!----------------------------------------------------------------------
9#if defined key_lim3
10   !!----------------------------------------------------------------------
11   !!   'key_lim3'                                       LIM3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!    ice_dia      : diagnostic of the sea-ice global heat content, salt content and volume conservation
14   !!    ice_dia_init : initialization of budget calculation
15   !!    ice_dia_rst  : read/write budgets restart
16   !!----------------------------------------------------------------------
17   USE ice            ! LIM-3: sea-ice variable
18   USE dom_oce        ! ocean domain
19   USE phycst         ! physical constant
20   USE daymod         ! model calendar
21   USE sbc_oce , ONLY : sfx, nn_fsbc   ! surface boundary condition: ocean fields
22   USE icerst         ! ice restart
23   !
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! MPP library
26   USE timing         ! preformance summary
27   USE iom            ! I/O manager
28   USE lib_fortran    ! glob_sum
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ice_dia        ! called by icestp.F90
34   PUBLIC   ice_dia_init   ! called in icestp.F90
35
36   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents
37   REAL(wp)                              ::   frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot  ! global forcing trends
38   
39   !! * Substitutions
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
43   !! $Id: icedia.F90 8413 2017-08-07 17:05:39Z clem $
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   INTEGER FUNCTION ice_dia_alloc()
49      !!---------------------------------------------------------------------!
50      !!                ***  ROUTINE ice_rdgrft_alloc ***
51      !!---------------------------------------------------------------------!
52      ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ice_dia_alloc )
53
54      IF( lk_mpp             )   CALL mpp_sum ( ice_dia_alloc )
55      IF( ice_dia_alloc /= 0 )   CALL ctl_warn( 'ice_dia_alloc: failed to allocate arrays' )
56      !
57   END FUNCTION ice_dia_alloc
58
59   SUBROUTINE ice_dia( kt )
60      !!---------------------------------------------------------------------------
61      !!                  ***  ROUTINE ice_dia  ***
62      !!     
63      !! ** Purpose:   Compute the sea-ice global heat content, salt content
64      !!             and volume conservation
65      !!---------------------------------------------------------------------------
66      INTEGER, INTENT(in) ::   kt   ! ocean time step
67      !!
68      REAL(wp)   ::   zbg_ivol, zbg_item, zbg_area, zbg_isal
69      REAL(wp)   ::   zbg_svol, zbg_stem
70      REAL(wp)   ::   z_frc_voltop, z_frc_temtop, z_frc_sal
71      REAL(wp)   ::   z_frc_volbot, z_frc_tembot 
72      REAL(wp)   ::   zdiff_vol, zdiff_sal, zdiff_tem 
73      !!---------------------------------------------------------------------------
74      IF( nn_timing == 1 )   CALL timing_start('ice_dia')
75
76      IF( kt == nit000 .AND. lwp ) THEN
77         WRITE(numout,*)
78         WRITE(numout,*)'icedia : outpout ice diagnostics (integrated over the domain)'
79         WRITE(numout,*)'~~~~~~'
80      ENDIF
81
82!!gm glob_sum includes a " * tmask_i ", so remove  " * tmask(:,:,1) "
83
84      ! ----------------------- !
85      ! 1 -  Contents !
86      ! ----------------------- !
87      zbg_ivol = glob_sum( vt_i(:,:) * e1e2t(:,:) ) * 1.e-9                  ! ice volume (km3)
88      zbg_svol = glob_sum( vt_s(:,:) * e1e2t(:,:) ) * 1.e-9                  ! snow volume (km3)
89      zbg_area = glob_sum( at_i(:,:) * e1e2t(:,:) ) * 1.e-6                  ! area (km2)
90      zbg_isal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e1e2t(:,:) ) * 1.e-9 ! salt content (pss*km3)
91      zbg_item = glob_sum( et_i * e1e2t(:,:) ) * 1.e-20                      ! heat content (1.e20 J)
92      zbg_stem = glob_sum( et_s * e1e2t(:,:) ) * 1.e-20                      ! heat content (1.e20 J)
93     
94      ! ---------------------------!
95      ! 2 - Trends due to forcing  !
96      ! ---------------------------!
97      z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9  ! freshwater flux ice/snow-ocean
98      z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ) * 1.e-9                     ! freshwater flux ice/snow-atm
99      z_frc_sal    = r1_rau0 * glob_sum(   - sfx(:,:) * e1e2t(:,:) ) * 1.e-9                                          ! salt fluxes ice/snow-ocean
100      z_frc_tembot =           glob_sum( hfx_out(:,:) * e1e2t(:,:) ) * 1.e-20                                         ! heat on top of ocean (and below ice)
101      z_frc_temtop =           glob_sum( hfx_in (:,:) * e1e2t(:,:) ) * 1.e-20                                         ! heat on top of ice-coean
102      !
103      frc_voltop  = frc_voltop  + z_frc_voltop  * rdt_ice ! km3
104      frc_volbot  = frc_volbot  + z_frc_volbot  * rdt_ice ! km3
105      frc_sal     = frc_sal     + z_frc_sal     * rdt_ice ! km3*pss
106      frc_temtop  = frc_temtop  + z_frc_temtop  * rdt_ice ! 1.e20 J
107      frc_tembot  = frc_tembot  + z_frc_tembot  * rdt_ice ! 1.e20 J
108           
109      ! ----------------------- !
110      ! 3 -  Content variations !
111      ! ----------------------- !
112      zdiff_vol = r1_rau0 * glob_sum( ( rhoic*vt_i(:,:) + rhosn*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3)
113      zdiff_sal = r1_rau0 * glob_sum( ( rhoic* SUM( smv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! salt content trend (km3*pss)
114      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:)             - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20 ! heat content trend (1.e20 J)
115      !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check)
116
117      ! ----------------------- !
118      ! 4 -  Drifts             !
119      ! ----------------------- !
120      zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot )
121      zdiff_sal = zdiff_sal - frc_sal
122      zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop )
123
124      ! ----------------------- !
125      ! 5 - Diagnostics writing !
126      ! ----------------------- !
127!!gm I don't understand the division by the ocean surface (i.e. glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt )
128!!   and its multiplication bu kt ! is it really what we want ? what is this quantity ?
129!!   IF it is really what we want, compute it at kt=nit000, not 3 time by time-step !
130!!   kt*rdt  : you mean rdtice ?
131!!gm
132      !
133      IF( iom_use('ibgvolume')    )   CALL iom_put( 'ibgvolume' , zdiff_vol     )   ! ice/snow volume  drift            (km3 equivalent ocean water)         
134      IF( iom_use('ibgsaltco')    )   CALL iom_put( 'ibgsaltco' , zdiff_sal     )   ! ice salt content drift            (psu*km3 equivalent ocean water)
135      IF( iom_use('ibgheatco')    )   CALL iom_put( 'ibgheatco' , zdiff_tem     )   ! ice/snow heat content drift       (1.e20 J)
136      IF( iom_use('ibgheatfx')    )   CALL iom_put( 'ibgheatfx' ,               &   ! ice/snow heat flux drift          (W/m2)
137         &                                                     zdiff_tem /glob_sum( e1e2t(:,:) * 1.e-20 * kt*rdt ) )
138
139      IF( iom_use('ibgfrcvoltop') )   CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water)
140      IF( iom_use('ibgfrcvolbot') )   CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water)
141      IF( iom_use('ibgfrcsal')    )   CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)   
142      IF( iom_use('ibgfrctemtop') )   CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)   
143      IF( iom_use('ibgfrctembot') )   CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)   
144      IF( iom_use('ibgfrchfxtop') )   CALL iom_put( 'ibgfrchfxtop' ,            &   ! heat on top of ice/snw/ocean      (W/m2)
145         &                                                          frc_temtop / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt  )
146      IF( iom_use('ibgfrchfxbot') )   CALL iom_put( 'ibgfrchfxbot' ,            &   ! heat on top of ocean(below ice)   (W/m2)
147         &                                                          frc_tembot / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt  )
148
149      IF( iom_use('ibgvol_tot' )  )   CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                       (km3)
150      IF( iom_use('sbgvol_tot' )  )   CALL iom_put( 'sbgvol_tot'  , zbg_svol     )   ! snow volume                      (km3)
151      IF( iom_use('ibgarea_tot')  )   CALL iom_put( 'ibgarea_tot' , zbg_area     )   ! ice area                         (km2)
152      IF( iom_use('ibgsalt_tot')  )   CALL iom_put( 'ibgsalt_tot' , zbg_isal     )   ! ice salinity content             (pss*km3)
153      IF( iom_use('ibgheat_tot')  )   CALL iom_put( 'ibgheat_tot' , zbg_item     )   ! ice heat content                 (1.e20 J)
154      IF( iom_use('sbgheat_tot')  )   CALL iom_put( 'sbgheat_tot' , zbg_stem     )   ! snow heat content                (1.e20 J)
155      !
156      IF( lrst_ice )   CALL ice_dia_rst( 'WRITE', kt_ice )
157      !
158      IF( nn_timing == 1 )   CALL timing_stop('ice_dia')
159      !
160   END SUBROUTINE ice_dia
161
162
163   SUBROUTINE ice_dia_init
164      !!---------------------------------------------------------------------------
165      !!                  ***  ROUTINE ice_dia_init  ***
166      !!     
167      !! ** Purpose: Initialization for the heat salt volume budgets
168      !!
169      !! ** Method : Compute initial heat content, salt content and volume
170      !!
171      !! ** Action : - Compute initial heat content, salt content and volume
172      !!             - Initialize forcing trends
173      !!             - Compute coefficients for conversion
174      !!---------------------------------------------------------------------------
175      INTEGER            ::   ios, ierror   ! local integer
176      !!
177      NAMELIST/namice_dia/ ln_icediachk, ln_icediahsb, ln_icectl, iiceprt, jiceprt 
178      !!----------------------------------------------------------------------
179      !
180      REWIND( numnam_ice_ref )      ! Namelist namice_dia in reference namelist : Parameters for ice
181      READ  ( numnam_ice_ref, namice_dia, IOSTAT = ios, ERR = 901)
182901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_dia in reference namelist', lwp )
183
184      REWIND( numnam_ice_cfg )      ! Namelist namice_dia in configuration namelist : Parameters for ice
185      READ  ( numnam_ice_cfg, namice_dia, IOSTAT = ios, ERR = 902 )
186902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_dia in configuration namelist', lwp )
187      IF(lwm) WRITE ( numoni, namice_dia )
188      !
189      IF(lwp) THEN                  ! control print
190         WRITE(numout,*)
191         WRITE(numout,*) 'ice_dia_init: ice diagnostics'
192         WRITE(numout,*) ' ~~~~~~~~~~~'
193         WRITE(numout,*) '   Namelist namice_dia : '
194         WRITE(numout,*) '      Diagnose online heat/mass/salt budget      ln_icediachk = ', ln_icediachk
195         WRITE(numout,*) '      Output          heat/mass/salt budget      ln_icediahsb = ', ln_icediahsb
196         WRITE(numout,*) '      control prints for a given grid point         ln_icectl = ', ln_icectl
197         WRITE(numout,*) '         chosen grid point position         (iiceprt,jiceprt) = (', iiceprt,',', jiceprt,')'
198      ENDIF
199      !     
200      IF( ln_icediahsb ) THEN
201         IF( ice_dia_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dia_init : unable to allocate arrays' )   ! allocate tke arrays
202         CALL ice_dia_rst( 'READ' )   ! read or initialize all required files
203      ENDIF
204      !
205   END SUBROUTINE ice_dia_init
206
207
208   SUBROUTINE ice_dia_rst( cdrw, kt )
209      !!---------------------------------------------------------------------
210      !!                   ***  ROUTINE icedia_rst  ***
211      !!                     
212      !! ** Purpose :   Read or write DIA file in restart file
213      !!
214      !! ** Method  :   use of IOM library
215      !!----------------------------------------------------------------------
216      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
217      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
218      !
219      INTEGER  ::   iter    ! local integer
220      REAL(wp) ::   ziter   ! local scalar
221      !!----------------------------------------------------------------------
222      !
223      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
224         IF( ln_rstart ) THEN                   !* Read the restart file
225            !
226            CALL iom_get( numrir, 'kt_ice' , ziter )
227            IF(lwp) WRITE(numout,*)
228            IF(lwp) WRITE(numout,*) 'ice_dia_rst read at time step = ', ziter
229            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
230            CALL iom_get( numrir, 'frc_voltop' , frc_voltop  )
231            CALL iom_get( numrir, 'frc_volbot' , frc_volbot  )
232            CALL iom_get( numrir, 'frc_temtop' , frc_temtop  )
233            CALL iom_get( numrir, 'frc_tembot' , frc_tembot  )
234            CALL iom_get( numrir, 'frc_sal'    , frc_sal     )
235            CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini )
236            CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini )
237            CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini )
238         ELSE
239            IF(lwp) WRITE(numout,*)
240            IF(lwp) WRITE(numout,*) ' ice_dia at initial state '
241            IF(lwp) WRITE(numout,*) '~~~~~~~'
242            ! set trends to 0
243            frc_voltop  = 0._wp                                         
244            frc_volbot  = 0._wp                                         
245            frc_temtop  = 0._wp                                                 
246            frc_tembot  = 0._wp                                                 
247            frc_sal     = 0._wp                                                 
248            ! record initial ice volume, salt and temp
249            vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2)
250            tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J)
251            sal_loc_ini(:,:) = rhoic * SUM( smv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2)
252         ENDIF
253         !
254      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
255         !                                   ! -------------------
256         iter = kt + nn_fsbc - 1   ! ice restarts are written at kt == nitrst - nn_fsbc + 1
257         !
258         IF( iter == nitrst ) THEN
259            IF(lwp) WRITE(numout,*)
260            IF(lwp) WRITE(numout,*) 'ice_dia_rst write at time step = ', kt
261            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
262         ENDIF
263         !
264         ! Write in numriw (if iter == nitrst)
265         ! ------------------
266         CALL iom_rstput( iter, nitrst, numriw, 'frc_voltop' , frc_voltop  )
267         CALL iom_rstput( iter, nitrst, numriw, 'frc_volbot' , frc_volbot  )
268         CALL iom_rstput( iter, nitrst, numriw, 'frc_temtop' , frc_temtop  )
269         CALL iom_rstput( iter, nitrst, numriw, 'frc_tembot' , frc_tembot  )
270         CALL iom_rstput( iter, nitrst, numriw, 'frc_sal'    , frc_sal     )
271         CALL iom_rstput( iter, nitrst, numriw, 'vol_loc_ini', vol_loc_ini )
272         CALL iom_rstput( iter, nitrst, numriw, 'tem_loc_ini', tem_loc_ini )
273         CALL iom_rstput( iter, nitrst, numriw, 'sal_loc_ini', sal_loc_ini )
274         !
275      ENDIF
276      !
277   END SUBROUTINE ice_dia_rst
278 
279#else
280   !!----------------------------------------------------------------------
281   !!   Default option :         Empty module          NO LIM sea-ice model
282   !!----------------------------------------------------------------------
283#endif
284
285   !!======================================================================
286END MODULE icedia
Note: See TracBrowser for help on using the repository browser.