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.
icedia.F90 in NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/ICE/icedia.F90 @ 11954

Last change on this file since 11954 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
  • Property svn:keywords set to Id
File size: 15.2 KB
Line 
1MODULE icedia
2   !!======================================================================
3   !!                       ***  MODULE icedia  ***
4   !!  Sea-Ice:   global budgets
5   !!======================================================================
6   !! History :  3.4  !  2012-10  (C. Rousset)       original code
7   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 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 dom_oce        ! ocean domain
18   USE phycst         ! physical constant
19   USE daymod         ! model calendar
20   USE sbc_oce , ONLY : sfx, nn_fsbc   ! surface boundary condition: ocean fields
21   USE ice            ! sea-ice: variables
22   USE icerst         ! sea-ice: restart
23   !
24   USE in_out_manager ! I/O manager
25   USE iom            ! I/O manager library
26   USE lib_mpp        ! MPP library
27   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
28   USE timing         ! Timing
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), SAVE ::   z1_e1e2  ! inverse of the ocean area
37   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini                    ! initial volume, salt and heat contents
38   REAL(wp)                              ::   frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot  ! global forcing trends
39   
40   !! * Substitutions
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
44   !! $Id$
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   INTEGER FUNCTION ice_dia_alloc()
50      !!---------------------------------------------------------------------!
51      !!                ***  ROUTINE ice_dia_alloc ***
52      !!---------------------------------------------------------------------!
53      ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ice_dia_alloc )
54
55      CALL mpp_sum ( 'icedia', ice_dia_alloc )
56      IF( ice_dia_alloc /= 0 )   CALL ctl_stop( 'STOP',  'ice_dia_alloc: failed to allocate arrays'  )
57      !
58   END FUNCTION ice_dia_alloc
59
60   SUBROUTINE ice_dia( kt )
61      !!---------------------------------------------------------------------------
62      !!                  ***  ROUTINE ice_dia  ***
63      !!     
64      !! ** Purpose:   Compute the sea-ice global heat content, salt content
65      !!             and volume conservation
66      !!---------------------------------------------------------------------------
67      INTEGER, INTENT(in) ::   kt   ! ocean time step
68      !!
69      REAL(wp)   ::   zbg_ivol, zbg_item, zbg_area, zbg_isal
70      REAL(wp)   ::   zbg_svol, zbg_stem
71      REAL(wp)   ::   z_frc_voltop, z_frc_temtop, z_frc_sal
72      REAL(wp)   ::   z_frc_volbot, z_frc_tembot 
73      REAL(wp)   ::   zdiff_vol, zdiff_sal, zdiff_tem 
74      !!---------------------------------------------------------------------------
75      IF( ln_timing )   CALL timing_start('ice_dia')
76
77      IF( kt == nit000 .AND. lwp ) THEN
78         WRITE(numout,*)
79         WRITE(numout,*)'icedia: output ice diagnostics (integrated over the domain)'
80         WRITE(numout,*)'~~~~~~'
81      ENDIF
82
83      IF( kt == nit000 ) THEN
84         z1_e1e2 = 1._wp / glob_sum( 'icedia', e1e2t(:,:) )
85      ENDIF
86     
87      ! ----------------------- !
88      ! 1 -  Contents           !
89      ! ----------------------- !
90      IF(  iom_use('ibgvol_tot' ) .OR. iom_use('sbgvol_tot' ) .OR. iom_use('ibgarea_tot') .OR. &
91         & iom_use('ibgsalt_tot') .OR. iom_use('ibgheat_tot') .OR. iom_use('sbgheat_tot') ) THEN
92
93         zbg_ivol = glob_sum( 'icedia', vt_i(:,:) * e1e2t(:,:) ) * 1.e-9  ! ice volume (km3)
94         zbg_svol = glob_sum( 'icedia', vt_s(:,:) * e1e2t(:,:) ) * 1.e-9  ! snow volume (km3)
95         zbg_area = glob_sum( 'icedia', at_i(:,:) * e1e2t(:,:) ) * 1.e-6  ! area (km2)
96         zbg_isal = glob_sum( 'icedia', st_i(:,:) * e1e2t(:,:) ) * 1.e-9  ! salt content (pss*km3)
97         zbg_item = glob_sum( 'icedia', et_i(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J)
98         zbg_stem = glob_sum( 'icedia', et_s(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J)
99
100         CALL iom_put( 'ibgvol_tot'  , zbg_ivol ) 
101         CALL iom_put( 'sbgvol_tot'  , zbg_svol ) 
102         CALL iom_put( 'ibgarea_tot' , zbg_area ) 
103         CALL iom_put( 'ibgsalt_tot' , zbg_isal ) 
104         CALL iom_put( 'ibgheat_tot' , zbg_item ) 
105         CALL iom_put( 'sbgheat_tot' , zbg_stem ) 
106 
107      ENDIF
108
109      ! ---------------------------!
110      ! 2 - Trends due to forcing  !
111      ! ---------------------------!
112      ! they must be kept outside an IF(iom_use) because of the call to dia_rst below
113      z_frc_volbot = r1_rau0 * glob_sum( 'icedia', -( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-ocean
114      z_frc_voltop = r1_rau0 * glob_sum( 'icedia', -( wfx_sub(:,:) + wfx_spr(:,:) )                    * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-atm
115      z_frc_sal    = r1_rau0 * glob_sum( 'icedia', -      sfx(:,:)                                     * e1e2t(:,:) ) * 1.e-9   ! salt fluxes ice/snow-ocean
116      z_frc_tembot =           glob_sum( 'icedia',  qt_oce_ai(:,:)                                     * e1e2t(:,:) ) * 1.e-20  ! heat on top of ocean (and below ice)
117      z_frc_temtop =           glob_sum( 'icedia',  qt_atm_oi(:,:)                                     * e1e2t(:,:) ) * 1.e-20  ! heat on top of ice-coean
118      !
119      frc_voltop  = frc_voltop  + z_frc_voltop  * rdt_ice ! km3
120      frc_volbot  = frc_volbot  + z_frc_volbot  * rdt_ice ! km3
121      frc_sal     = frc_sal     + z_frc_sal     * rdt_ice ! km3*pss
122      frc_temtop  = frc_temtop  + z_frc_temtop  * rdt_ice ! 1.e20 J
123      frc_tembot  = frc_tembot  + z_frc_tembot  * rdt_ice ! 1.e20 J
124
125      CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water)
126      CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water)
127      CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)   
128      CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)   
129      CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)   
130
131      IF(  iom_use('ibgfrchfxtop') .OR. iom_use('ibgfrchfxbot') ) THEN
132         CALL iom_put( 'ibgfrchfxtop' , frc_temtop * z1_e1e2 * 1.e-20 * kt*rdt ) ! heat on top of ice/snw/ocean      (W/m2)
133         CALL iom_put( 'ibgfrchfxbot' , frc_tembot * z1_e1e2 * 1.e-20 * kt*rdt ) ! heat on top of ocean(below ice)   (W/m2)
134      ENDIF
135     
136      ! ---------------------------------- !
137      ! 3 -  Content variations and drifts !
138      ! ---------------------------------- !
139      IF(  iom_use('ibgvolume') .OR. iom_use('ibgsaltco') .OR. iom_use('ibgheatco') .OR. iom_use('ibgheatfx') ) THEN
140           
141         zdiff_vol = r1_rau0 * glob_sum( 'icedia', ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater trend (km3)
142         zdiff_sal = r1_rau0 * glob_sum( 'icedia', ( rhoi*st_i(:,:)                  - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! salt content trend (km3*pss)
143         zdiff_tem =           glob_sum( 'icedia', ( et_i(:,:) + et_s(:,:)           - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20  ! heat content trend (1.e20 J)
144         !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check)
145         
146         zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot )
147         zdiff_sal = zdiff_sal - frc_sal
148         zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop )
149         
150         CALL iom_put( 'ibgvolume' , zdiff_vol )   ! ice/snow volume  drift            (km3 equivalent ocean water)         
151         CALL iom_put( 'ibgsaltco' , zdiff_sal )   ! ice salt content drift            (psu*km3 equivalent ocean water)
152         CALL iom_put( 'ibgheatco' , zdiff_tem )   ! ice/snow heat content drift       (1.e20 J)
153         !
154      ENDIF
155     
156      IF( lrst_ice )   CALL ice_dia_rst( 'WRITE', kt_ice )
157      !
158      IF( ln_timing )   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/namdia/ ln_icediachk, rn_icechk_cel, rn_icechk_glo, ln_icediahsb, ln_icectl, iiceprt, jiceprt 
178      !!----------------------------------------------------------------------
179      !
180      READ  ( numnam_ice_ref, namdia, IOSTAT = ios, ERR = 901)
181901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdia in reference namelist' )
182      READ  ( numnam_ice_cfg, namdia, IOSTAT = ios, ERR = 902 )
183902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdia in configuration namelist' )
184      IF(lwm) WRITE ( numoni, namdia )
185      !
186      IF(lwp) THEN                  ! control print
187         WRITE(numout,*)
188         WRITE(numout,*) 'ice_dia_init: ice diagnostics'
189         WRITE(numout,*) ' ~~~~~~~~~~~'
190         WRITE(numout,*) '   Namelist namdia:'
191         WRITE(numout,*) '      Diagnose online heat/mass/salt conservation ln_icediachk  = ', ln_icediachk
192         WRITE(numout,*) '         threshold for conservation (gridcell)    rn_icechk_cel = ', rn_icechk_cel
193         WRITE(numout,*) '         threshold for conservation (global)      rn_icechk_glo = ', rn_icechk_glo
194         WRITE(numout,*) '      Output          heat/mass/salt budget       ln_icediahsb  = ', ln_icediahsb
195         WRITE(numout,*) '      control prints for a given grid point       ln_icectl     = ', ln_icectl
196         WRITE(numout,*) '         chosen grid point position          (iiceprt,jiceprt)  = (', iiceprt,',', jiceprt,')'
197      ENDIF
198      !     
199      IF( ln_icediahsb ) THEN
200         IF( ice_dia_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dia_init : unable to allocate arrays' )   ! allocate tke arrays
201         CALL ice_dia_rst( 'READ' )   ! read or initialize all required files
202      ENDIF
203      !
204   END SUBROUTINE ice_dia_init
205
206
207   SUBROUTINE ice_dia_rst( cdrw, kt )
208      !!---------------------------------------------------------------------
209      !!                   ***  ROUTINE icedia_rst  ***
210      !!                     
211      !! ** Purpose :   Read or write DIA file in restart file
212      !!
213      !! ** Method  :   use of IOM library
214      !!----------------------------------------------------------------------
215      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
216      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
217      !
218      INTEGER  ::   iter    ! local integer
219      REAL(wp) ::   ziter   ! local scalar
220      !!----------------------------------------------------------------------
221      !
222      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
223         IF( ln_rstart ) THEN                   !* Read the restart file
224            !
225            CALL iom_get( numrir, 'kt_ice' , ziter )
226            IF(lwp) WRITE(numout,*)
227            IF(lwp) WRITE(numout,*) 'ice_dia_rst read at time step = ', ziter
228            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
229            CALL iom_get( numrir, 'frc_voltop' , frc_voltop  )
230            CALL iom_get( numrir, 'frc_volbot' , frc_volbot  )
231            CALL iom_get( numrir, 'frc_temtop' , frc_temtop  )
232            CALL iom_get( numrir, 'frc_tembot' , frc_tembot  )
233            CALL iom_get( numrir, 'frc_sal'    , frc_sal     )
234            CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini )
235            CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini )
236            CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini )
237         ELSE
238            IF(lwp) WRITE(numout,*)
239            IF(lwp) WRITE(numout,*) ' ice_dia at initial state '
240            IF(lwp) WRITE(numout,*) '~~~~~~~'
241            ! set trends to 0
242            frc_voltop  = 0._wp                                         
243            frc_volbot  = 0._wp                                         
244            frc_temtop  = 0._wp                                                 
245            frc_tembot  = 0._wp                                                 
246            frc_sal     = 0._wp                                                 
247            ! record initial ice volume, salt and temp
248            vol_loc_ini(:,:) = rhoi * vt_i(:,:) + rhos * vt_s(:,:)  ! ice/snow volume (kg/m2)
249            tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                ! ice/snow heat content (J)
250            sal_loc_ini(:,:) = rhoi * st_i(:,:)                     ! ice salt content (pss*kg/m2)
251         ENDIF
252         !
253      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
254         !                                   ! -------------------
255         iter = kt + nn_fsbc - 1   ! ice restarts are written at kt == nitrst - nn_fsbc + 1
256         !
257         IF( iter == nitrst ) THEN
258            IF(lwp) WRITE(numout,*)
259            IF(lwp) WRITE(numout,*) 'ice_dia_rst write at time step = ', kt
260            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
261         ENDIF
262         !
263         ! Write in numriw (if iter == nitrst)
264         ! ------------------
265         CALL iom_rstput( iter, nitrst, numriw, 'frc_voltop' , frc_voltop  )
266         CALL iom_rstput( iter, nitrst, numriw, 'frc_volbot' , frc_volbot  )
267         CALL iom_rstput( iter, nitrst, numriw, 'frc_temtop' , frc_temtop  )
268         CALL iom_rstput( iter, nitrst, numriw, 'frc_tembot' , frc_tembot  )
269         CALL iom_rstput( iter, nitrst, numriw, 'frc_sal'    , frc_sal     )
270         CALL iom_rstput( iter, nitrst, numriw, 'vol_loc_ini', vol_loc_ini )
271         CALL iom_rstput( iter, nitrst, numriw, 'tem_loc_ini', tem_loc_ini )
272         CALL iom_rstput( iter, nitrst, numriw, 'sal_loc_ini', sal_loc_ini )
273         !
274      ENDIF
275      !
276   END SUBROUTINE ice_dia_rst
277 
278#else
279   !!----------------------------------------------------------------------
280   !!   Default option :         Empty module         NO SI3 sea-ice model
281   !!----------------------------------------------------------------------
282#endif
283
284   !!======================================================================
285END MODULE icedia
Note: See TracBrowser for help on using the repository browser.