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 NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/tests/STATION_ASF/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/tests/STATION_ASF/MY_SRC/icestp.F90 @ 13655

Last change on this file since 13655 was 13655, checked in by laurent, 4 years ago

Commit all my dev of 2020!
File size: 16.8 KB
Line 
1MODULE icestp
2   !!======================================================================
3   !!                       ***  MODULE  icestp  ***
4   !! sea ice : Master routine for all the sea ice model
5   !!=====================================================================
6   !!
7   !! The sea ice model SI3 (Sea Ice modelling Integrated Initiative),
8   !!                        aka Sea Ice cube for its nickname
9   !!
10   !!    is originally based on LIM3, developed in Louvain-la-Neuve by:
11   !!       * Martin Vancoppenolle (UCL-ASTR, Belgium)
12   !!       * Sylvain Bouillon (UCL-ASTR, Belgium)
13   !!       * Miguel Angel Morales Maqueda (NOC-L, UK)
14   !!      thanks to valuable earlier work by
15   !!       * Thierry Fichefet
16   !!       * Hugues Goosse
17   !!      thanks also to the following persons who contributed
18   !!       * Gurvan Madec, Claude Talandier, Christian Ethe (LOCEAN, France)
19   !!       * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany)
20   !!       * Bill Lipscomb (LANL), Cecilia Bitz (UWa) and Elisabeth Hunke (LANL), USA.
21   !!
22   !! SI3 has been made possible by a handful of persons who met as working group
23   !!      (from France, Belgium, UK and Italy)
24   !!    * Clement Rousset, Martin Vancoppenolle & Gurvan Madec (LOCEAN, France)
25   !!    * Matthieu Chevalier & David Salas (Meteo France, France)
26   !!    * Gilles Garric (Mercator Ocean, France)
27   !!    * Thierry Fichefet & Francois Massonnet (UCL, Belgium)
28   !!    * Ed Blockley & Jeff Ridley (Met Office, UK)
29   !!    * Danny Feltham & David Schroeder (CPOM, UK)
30   !!    * Yevgeny Aksenov (NOC, UK)
31   !!    * Paul Holland (BAS, UK)
32   !!    * Dorotea Iovino (CMCC, Italy)
33   !!======================================================================
34   !! History :  4.0  !  2018     (C. Rousset)      Original code SI3
35   !!----------------------------------------------------------------------
36#if defined key_si3
37   !!----------------------------------------------------------------------
38   !!   'key_si3'                                       SI3 sea-ice model
39   !!----------------------------------------------------------------------
40   !!   ice_stp       : sea-ice model time-stepping and update ocean SBC over ice-covered area
41   !!   ice_init      : initialize sea-ice
42   !!----------------------------------------------------------------------
43   USE oce            ! ocean dynamics and tracers
44   USE dom_oce        ! ocean space and time domain
45   USE c1d            ! 1D vertical configuration
46   USE ice            ! sea-ice: variables
47   USE ice1D          ! sea-ice: thermodynamical 1D variables
48   !
49   USE phycst         ! Define parameters for the routines
50   USE eosbn2         ! equation of state
51   USE sbc_oce        ! Surface boundary condition: ocean fields
52   USE sbc_ice        ! Surface boundary condition: ice   fields
53   !
54   USE icesbc         ! sea-ice: Surface boundary conditions
55   USE iceupdate      ! sea-ice: sea surface boundary condition update
56   USE icewri         ! sea-ice: outputs
57   !
58   !
59   USE in_out_manager ! I/O manager
60   USE iom            ! I/O manager library
61   USE lib_mpp        ! MPP library
62   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
63   USE timing         ! Timing
64   USE prtctl         ! Print control
65
66   IMPLICIT NONE
67   PRIVATE
68
69   PUBLIC   ice_stp    ! called by sbcmod.F90
70   PUBLIC   ice_init   ! called by sbcmod.F90
71
72   !! * Substitutions
73#  include "do_loop_substitute.h90"
74   !!----------------------------------------------------------------------
75   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
76   !! $Id: icestp.F90 13641 2020-10-19 18:16:58Z clem $
77   !! Software governed by the CeCILL license (see ./LICENSE)
78   !!----------------------------------------------------------------------
79CONTAINS
80
81   SUBROUTINE ice_stp( kt, Kbb, Kmm, ksbc )
82      !!---------------------------------------------------------------------
83      !!                  ***  ROUTINE ice_stp  ***
84      !!
85      !! ** Purpose :   sea-ice model time-stepping and update ocean surface
86      !!              boundary condition over ice-covered area
87      !!
88      !! ** Method  :   ice model time stepping
89      !!              - call the ice dynamics routine
90      !!              - call the ice advection/diffusion routine
91      !!              - call the ice thermodynamics routine
92      !!              - call the routine that computes mass and
93      !!                heat fluxes at the ice/ocean interface
94      !!              - save the outputs
95      !!              - save the outputs for restart when necessary
96      !!
97      !! ** Action  : - time evolution of the LIM sea-ice model
98      !!              - update all sbc variables below sea-ice:
99      !!                utau, vtau, taum, wndm, qns , qsr, emp , sfx
100      !!---------------------------------------------------------------------
101      INTEGER, INTENT(in) ::   kt       ! ocean time step
102      INTEGER, INTENT(in) ::   Kbb, Kmm ! ocean time level indices
103      INTEGER, INTENT(in) ::   ksbc     ! flux formulation (user defined, bulk, or Pure Coupled)
104      !
105      INTEGER ::   jl   ! dummy loop index
106      !!----------------------------------------------------------------------
107      !
108      IF( ln_timing )   CALL timing_start('ice_stp')
109      !
110      !                                      !-----------------------!
111      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN   ! --- Ice time step --- !
112         !                                   !-----------------------!
113         !
114         kt_ice = kt                              ! -- Ice model time step
115         !
116         u_oce(:,:) = ssu_m(:,:)                  ! -- mean surface ocean current
117         v_oce(:,:) = ssv_m(:,:)
118         !
119         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )   ! -- freezing temperature [Kelvin] (set to rt0 over land)
120         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
121         !
122         !
123         !------------------------------------------------!
124         ! --- Dynamical coupling with the atmosphere --- !
125         !------------------------------------------------!
126         ! It provides the following fields used in sea ice model:
127         !    utau_ice, vtau_ice = surface ice stress [N/m2]
128         !------------------------------------------------!
129         !------------------------------------------------------!
130         ! --- Thermodynamical coupling with the atmosphere --- !
131         !------------------------------------------------------!
132         ! It provides the following fields used in sea ice model:
133         !    emp_oce , emp_ice    = E-P over ocean and sea ice                    [Kg/m2/s]
134         !    sprecip              = solid precipitation                           [Kg/m2/s]
135         !    evap_ice             = sublimation                                   [Kg/m2/s]
136         !    qsr_tot , qns_tot    = solar & non solar heat flux (total)           [W/m2]
137         !    qsr_ice , qns_ice    = solar & non solar heat flux over ice          [W/m2]
138         !    dqns_ice             = non solar  heat sensistivity                  [W/m2]
139         !    qemp_oce, qemp_ice,  = sensible heat (associated with evap & precip) [W/m2]
140         !    qprec_ice, qevap_ice
141         !------------------------------------------------------!
142                                        CALL ice_sbc_flx( kt, ksbc )
143         !----------------------------!
144         ! --- ice thermodynamics --- !
145         !----------------------------!
146
147         fr_i  (:,:)   = at_i(:,:) !#LB...
148
149         !!#LB: lines stolen from "ice_update_flx()@iceupdate.F90"
150         !!=============================================================
151
152         !IF(lwp) WRITE(numout,*) ''
153         !IF(lwp) WRITE(numout,*) 'LOLO:ice_stp()@icestp.F90 => calling "ice_update_flx()" at kt =', kt
154         !                               CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes
155         !IF(lwp) WRITE(numout,*) 'LOLO:ice_stp()@icestp.F90 => Done with "ice_update_flx()" at kt =', kt
156         !IF(lwp) WRITE(numout,*) ''
157
158         ! --- case we bypass ice thermodynamics --- !
159         !IF( .NOT. ln_icethd ) THEN   ! we suppose ice is impermeable => ocean is isolated from atmosphere
160         qt_atm_oi  (:,:)   = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
161         qt_oce_ai  (:,:)   = ( 1._wp - at_i_b(:,:) ) *   qns_oce(:,:)                  + qemp_oce(:,:)
162         emp_ice    (:,:)   = 0._wp
163         qemp_ice   (:,:)   = 0._wp
164         qevap_ice  (:,:,:) = 0._wp
165         !ENDIF
166
167         ! output all fluxes
168         !------------------
169
170         ! --- mass fluxes [kg/m2/s] --- !
171         IF( iom_use('emp_oce'    ) )   CALL iom_put( 'emp_oce', emp_oce )   ! emp over ocean (taking into account the snow blown away from the ice)
172         IF( iom_use('emp_ice'    ) )   CALL iom_put( 'emp_ice', emp_ice )   ! emp over ice   (taking into account the snow blown away from the ice)
173         ! --- heat fluxes [W/m2] --- !
174         IF( iom_use('qsr_oce'    ) )   CALL iom_put( 'qsr_oce'    , qsr_oce * ( 1._wp - at_i_b )                               )   !     solar flux at ocean surface
175         IF( iom_use('qns_oce'    ) )   CALL iom_put( 'qns_oce'    , qns_oce * ( 1._wp - at_i_b ) + qemp_oce                    )   ! non-solar flux at ocean surface
176         IF( iom_use('qns_ice'    ) )   CALL iom_put( 'qns_ice'    , SUM( qns_ice * a_i_b, dim=3 ) + qemp_ice                   )   ! non-solar flux at ice surface
177         IF( iom_use('qsr_ice'    ) )   CALL iom_put( 'qsr_ice'    , SUM( qsr_ice * a_i_b, dim=3 )                              )
178         IF( iom_use('qt_oce'     ) )   CALL iom_put( 'qt_oce'     ,      ( qsr_oce + qns_oce ) * ( 1._wp - at_i_b ) + qemp_oce )
179         IF( iom_use('qt_ice'     ) )   CALL iom_put( 'qt_ice'     , SUM( ( qns_ice + qsr_ice ) * a_i_b, dim=3 )     + qemp_ice )
180         IF( iom_use('qemp_oce'   ) )   CALL iom_put( 'qemp_oce'   , qemp_oce                                                   )   ! Downward Heat Flux from E-P over ocean
181         IF( iom_use('qemp_ice'   ) )   CALL iom_put( 'qemp_ice'   , qemp_ice                                                   )   ! Downward Heat Flux from E-P over ice
182
183
184         !!=============================================================
185
186         !
187                                        CALL ice_wri( kt )            ! -- Ice outputs
188         !
189         !IF( lrst_ice )                 CALL ice_rst_write( kt )      ! -- Ice restart file
190         !
191         !IF( ln_icectl )                CALL ice_ctl( kt )            ! -- Control checks
192         !
193      ENDIF   ! End sea-ice time step only
194
195      !-------------------------!
196      ! --- Ocean time step --- !
197      !-------------------------!
198      !IF( ln_icedyn )                   CALL ice_update_tau( kt, uu(:,:,1,Kbb), vv(:,:,1,Kbb) )   ! -- update surface ocean stresses
199      !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
200      !
201      IF( ln_timing )   CALL timing_stop('ice_stp')
202      !
203   END SUBROUTINE ice_stp
204
205
206   SUBROUTINE ice_init( Kbb, Kmm, Kaa )
207      !!----------------------------------------------------------------------
208      !!                  ***  ROUTINE ice_init  ***
209      !!
210      !! ** purpose :   Initialize sea-ice parameters
211      !!----------------------------------------------------------------------
212      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
213      !
214      INTEGER ::   ierr
215      !!----------------------------------------------------------------------
216      IF(lwp) WRITE(numout,*)
217      IF(lwp) WRITE(numout,*) 'Sea Ice Model: SI3 (Sea Ice modelling Integrated Initiative)'
218      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
219      IF(lwp) WRITE(numout,*)
220      IF(lwp) WRITE(numout,*) 'ice_init: Arrays allocation & Initialization of all routines & init state'
221      IF(lwp) WRITE(numout,*) '~~~~~~~~'
222      !
223      !                                ! Load the reference and configuration namelist files and open namelist output file
224      CALL load_nml( numnam_ice_ref, 'namelist_ice_ref',    numout, lwm )
225      CALL load_nml( numnam_ice_cfg, 'namelist_ice_cfg',    numout, lwm )
226      IF(lwm) CALL ctl_opn( numoni , 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
227      !
228      CALL par_init                ! set some ice run parameters
229      !
230      !                                ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
231      ierr =        ice_alloc        ()      ! ice variables
232      ierr = ierr + sbc_ice_alloc    ()      ! surface boundary conditions
233      ierr = ierr + ice1D_alloc      ()      ! thermodynamics
234      !
235      CALL mpp_sum( 'icestp', ierr )
236      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
237      !
238      !
239      CALL ice_sbc_init                ! set ice-ocean and ice-atm. coupling parameters
240      !
241      IF( ln_rstart )   CALL iom_close( numrir )  ! close input ice restart file
242      !
243   END SUBROUTINE ice_init
244
245
246   SUBROUTINE par_init
247      !!-------------------------------------------------------------------
248      !!                  ***  ROUTINE par_init ***
249      !!
250      !! ** Purpose :   Definition generic parameters for ice model
251      !!
252      !! ** Method  :   Read namelist and check the parameter
253      !!                values called at the first timestep (nit000)
254      !!
255      !! ** input   :   Namelist nampar
256      !!-------------------------------------------------------------------
257      INTEGER  ::   ios                 ! Local integer
258      !!
259      NAMELIST/nampar/ jpl, nlay_i, nlay_s, ln_virtual_itd, ln_icedyn, ln_icethd, rn_amax_n, rn_amax_s,  &
260         &             cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir
261      !!-------------------------------------------------------------------
262      !
263      READ  ( numnam_ice_ref, nampar, IOSTAT = ios, ERR = 901)
264901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampar in reference namelist' )
265      READ  ( numnam_ice_cfg, nampar, IOSTAT = ios, ERR = 902 )
266902   IF( ios > 0 )   CALL ctl_nam ( ios , 'nampar in configuration namelist' )
267      IF(lwm) WRITE( numoni, nampar )
268      !
269      IF(lwp) THEN                  ! control print
270         WRITE(numout,*)
271         WRITE(numout,*) '   par_init: ice parameters shared among all the routines'
272         WRITE(numout,*) '   ~~~~~~~~'
273         WRITE(numout,*) '      Namelist nampar: '
274         WRITE(numout,*) '         number of ice  categories                           jpl       = ', jpl
275         WRITE(numout,*) '         number of ice  layers                               nlay_i    = ', nlay_i
276         WRITE(numout,*) '         number of snow layers                               nlay_s    = ', nlay_s
277         WRITE(numout,*) '         virtual ITD param for jpl=1 (T) or not (F)     ln_virtual_itd = ', ln_virtual_itd
278         WRITE(numout,*) '         Ice dynamics       (T) or not (F)                   ln_icedyn = ', ln_icedyn
279         WRITE(numout,*) '         Ice thermodynamics (T) or not (F)                   ln_icethd = ', ln_icethd
280         WRITE(numout,*) '         maximum ice concentration for NH                              = ', rn_amax_n
281         WRITE(numout,*) '         maximum ice concentration for SH                              = ', rn_amax_s
282      ENDIF
283      !                                        !--- change max ice concentration for roundoff errors
284      rn_amax_n = MIN( rn_amax_n, 1._wp - epsi10 )
285      rn_amax_s = MIN( rn_amax_s, 1._wp - epsi10 )
286      !                                        !--- check consistency
287      IF ( jpl > 1 .AND. ln_virtual_itd ) THEN
288         ln_virtual_itd = .FALSE.
289         IF(lwp) WRITE(numout,*)
290         IF(lwp) WRITE(numout,*) '   ln_virtual_itd forced to false as jpl>1, no need with multiple categories to emulate them'
291      ENDIF
292      !
293      IF( ln_cpl .AND. nn_cats_cpl /= 1 .AND. nn_cats_cpl /= jpl ) THEN
294         CALL ctl_stop( 'STOP', 'par_init: in coupled mode, nn_cats_cpl should be either 1 or jpl' )
295      ENDIF
296      !
297      rDt_ice   = REAL(nn_fsbc) * rn_Dt          !--- sea-ice timestep and its inverse
298      r1_Dt_ice = 1._wp / rDt_ice
299      IF(lwp) WRITE(numout,*)
300      IF(lwp) WRITE(numout,*) '      ice timestep rDt_ice = nn_fsbc*rn_Dt = ', rDt_ice
301      !
302      r1_nlay_i = 1._wp / REAL( nlay_i, wp )   !--- inverse of nlay_i and nlay_s
303      r1_nlay_s = 1._wp / REAL( nlay_s, wp )
304      !
305   END SUBROUTINE par_init
306
307#else
308   !!----------------------------------------------------------------------
309   !!   Default option           Dummy module         NO SI3 sea-ice model
310   !!----------------------------------------------------------------------
311CONTAINS
312   SUBROUTINE ice_stp ( kt, ksbc )     ! Dummy routine
313      INTEGER, INTENT(in) ::   kt, ksbc
314      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
315   END SUBROUTINE ice_stp
316   SUBROUTINE ice_init                 ! Dummy routine
317   END SUBROUTINE ice_init
318#endif
319
320   !!======================================================================
321END MODULE icestp
Note: See TracBrowser for help on using the repository browser.