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.
icesbc.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/icesbc.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: 17.1 KB
Line 
1MODULE icesbc
2   !!======================================================================
3   !!                       ***  MODULE  icesbc  ***
4   !! Sea-Ice :   air-ice sbc fields
5   !!=====================================================================
6   !! History :  4.0  !  2017-08  (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   USE oce            ! ocean dynamics and tracers
14   USE dom_oce        ! ocean space and time domain
15   USE ice            ! sea-ice: variables
16   USE sbc_oce        ! Surface boundary condition: ocean fields
17   USE sbc_ice        ! Surface boundary condition: ice   fields
18   USE usrdef_sbc     ! Surface boundary condition: user defined
19   USE sbcblk         ! Surface boundary condition: bulk
20   USE sbccpl         ! Surface boundary condition: coupled interface
21   USE icealb         ! sea-ice: albedo
22   !
23   USE in_out_manager ! I/O manager
24   USE iom            ! I/O manager library
25   USE lib_mpp        ! MPP library
26   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
27   USE lbclnk         ! lateral boundary conditions (or mpp links)
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC ice_sbc_tau   ! called by icestp.F90
34   PUBLIC ice_sbc_flx   ! called by icestp.F90
35   PUBLIC ice_sbc_init  ! called by icestp.F90
36
37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
41   !! $Id$
42   !! Software governed by the CeCILL license (see ./LICENSE)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice )
47      !!-------------------------------------------------------------------
48      !!                  ***  ROUTINE ice_sbc_tau  ***
49      !!
50      !! ** Purpose : provide surface boundary condition for sea ice (momentum)
51      !!
52      !! ** Action  : It provides the following fields:
53      !!              utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2]
54      !!-------------------------------------------------------------------
55      INTEGER                     , INTENT(in   ) ::   kt                   ! ocean time step
56      INTEGER                     , INTENT(in   ) ::   ksbc                 ! type of sbc flux
57      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   utau_ice, vtau_ice   ! air-ice stress   [N/m2]
58      !!
59      INTEGER  ::   ji, jj                 ! dummy loop index
60      REAL(wp), DIMENSION(jpi,jpj) ::   zutau_ice, zvtau_ice 
61      !!-------------------------------------------------------------------
62      !
63      IF( ln_timing )   CALL timing_start('ice_sbc')
64      !
65      IF( kt == nit000 .AND. lwp ) THEN
66         WRITE(numout,*)
67         WRITE(numout,*)'ice_sbc_tau: Surface boundary condition for sea ice (momentum)'
68         WRITE(numout,*)'~~~~~~~~~~~~~~~'
69      ENDIF
70      !
71      SELECT CASE( ksbc )
72         CASE( jp_usr     )   ;    CALL usrdef_sbc_ice_tau( kt )                 ! user defined formulation
73         CASE( jp_blk     )   ;    CALL blk_ice_tau                              ! Bulk         formulation
74         CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled      formulation
75      END SELECT
76      !
77      IF( ln_mixcpl) THEN                                                        ! Case of a mixed Bulk/Coupled formulation
78                                   CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
79         DO jj = 2, jpjm1
80            DO ji = 2, jpim1
81               utau_ice(ji,jj) = utau_ice(ji,jj) * xcplmask(ji,jj,0) + zutau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
82               vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
83            END DO
84         END DO
85         CALL lbc_lnk_multi( 'icesbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
86      ENDIF
87      !
88      IF( ln_timing )   CALL timing_stop('ice_sbc')
89      !
90   END SUBROUTINE ice_sbc_tau
91
92   
93   SUBROUTINE ice_sbc_flx( kt, ksbc )
94      !!-------------------------------------------------------------------
95      !!                  ***  ROUTINE ice_sbc_flx  ***
96      !!
97      !! ** Purpose : provide surface boundary condition for sea ice (flux)
98      !!
99      !! ** Action  : It provides the following fields used in sea ice model:
100      !!                emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s]
101      !!                sprecip                                  = solid precipitation                           [Kg/m2/s]
102      !!                evap_ice                                 = sublimation                                   [Kg/m2/s]
103      !!                qsr_tot , qns_tot                        = solar & non solar heat flux (total)           [W/m2]
104      !!                qsr_ice , qns_ice                        = solar & non solar heat flux over ice          [W/m2]
105      !!                dqns_ice                                 = non solar  heat sensistivity                  [W/m2]
106      !!                qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2]
107      !!            + some fields that are not used outside this module:
108      !!                qla_ice                                  = latent heat flux over ice                     [W/m2]
109      !!                dqla_ice                                 = latent heat sensistivity                      [W/m2]
110      !!                tprecip                                  = total  precipitation                          [Kg/m2/s]
111      !!                alb_ice                                  = albedo above sea ice
112      !!-------------------------------------------------------------------
113      INTEGER, INTENT(in) ::   kt     ! ocean time step
114      INTEGER, INTENT(in) ::   ksbc   ! flux formulation (user defined, bulk or Pure Coupled)
115      !
116      INTEGER  ::   ji, jj, jl      ! dummy loop index
117      REAL(wp) ::   zmiss_val       ! missing value retrieved from xios
118      REAL(wp), DIMENSION(jpi,jpj,jpl)              ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
119      REAL(wp), DIMENSION(:,:)        , ALLOCATABLE ::   zalb, zmsk00      ! 2D workspace
120      !!--------------------------------------------------------------------
121      !
122      IF( ln_timing )   CALL timing_start('ice_sbc_flx')
123
124      IF( kt == nit000 .AND. lwp ) THEN
125         WRITE(numout,*)
126         WRITE(numout,*)'ice_sbc_flx: Surface boundary condition for sea ice (flux)'
127         WRITE(numout,*)'~~~~~~~~~~~~~~~'
128      ENDIF
129
130      ! get missing value from xml
131      CALL iom_miss_val( "icetemp", zmiss_val )
132
133      ! --- cloud-sky and overcast-sky ice albedos --- !
134      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )
135
136      ! albedo depends on cloud fraction because of non-linear spectral effects
137!!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument !
138      alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
139      !
140      !
141      SELECT CASE( ksbc )   !== fluxes over sea ice ==!
142      !
143      CASE( jp_usr )              !--- user defined formulation
144                                  CALL usrdef_sbc_ice_flx( kt, h_s, h_i )
145      CASE( jp_blk )              !--- bulk formulation
146                                  CALL blk_ice_flx    ( t_su, h_s, h_i, alb_ice )    !
147         IF( ln_mixcpl        )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
148         IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
149         !                        !    compute conduction flux and surface temperature (as in Jules surface module)
150         IF( ln_cndflx .AND. .NOT.ln_cndemulate ) &
151            &                     CALL blk_ice_qcn    ( ln_virtual_itd, t_su, t_bo, h_s, h_i )
152      CASE ( jp_purecpl )         !--- coupled formulation
153                                  CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
154         IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
155      END SELECT
156
157      !--- output ice albedo and surface albedo ---!
158      IF( iom_use('icealb') .OR. iom_use('albedo') ) THEN
159
160         ALLOCATE( zalb(jpi,jpj), zmsk00(jpi,jpj) )
161
162         WHERE( at_i_b < 1.e-03 )
163            zmsk00(:,:) = 0._wp
164            zalb  (:,:) = rn_alb_oce
165         ELSEWHERE
166            zmsk00(:,:) = 1._wp           
167            zalb  (:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
168         END WHERE
169         ! ice albedo
170         CALL iom_put( 'icealb' , zalb * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )
171         ! ice+ocean albedo
172         zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b )
173         CALL iom_put( 'albedo' , zalb )
174
175         DEALLOCATE( zalb, zmsk00 )
176
177      ENDIF
178      !
179      IF( ln_timing )   CALL timing_stop('ice_sbc_flx')
180      !
181   END SUBROUTINE ice_sbc_flx
182
183
184   SUBROUTINE ice_flx_dist( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_flxdist )
185      !!-------------------------------------------------------------------
186      !!                  ***  ROUTINE ice_flx_dist  ***
187      !!
188      !! ** Purpose :   update the ice surface boundary condition by averaging
189      !!              and/or redistributing fluxes on ice categories
190      !!
191      !! ** Method  :   average then redistribute
192      !!
193      !! ** Action  :   depends on k_flxdist
194      !!                = -1  Do nothing (needs N(cat) fluxes)
195      !!                =  0  Average N(cat) fluxes then apply the average over the N(cat) ice
196      !!                =  1  Average N(cat) fluxes then redistribute over the N(cat) ice
197      !!                                                 using T-ice and albedo sensitivity
198      !!                =  2  Redistribute a single flux over categories
199      !!-------------------------------------------------------------------
200      INTEGER                   , INTENT(in   ) ::   k_flxdist  ! redistributor
201      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature
202      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo
203      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqns_ice   ! non solar flux
204      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux
205      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdqn_ice   ! non solar flux sensitivity
206      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pevap_ice  ! sublimation
207      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity
208      !
209      INTEGER  ::   jl      ! dummy loop index
210      !
211      REAL(wp), DIMENSION(jpi,jpj) ::   z1_at_i   ! inverse of concentration
212      !
213      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qsr_m   ! Mean solar heat flux over all categories
214      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qns_m   ! Mean non solar heat flux over all categories
215      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_evap_m  ! Mean sublimation over all categories
216      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_dqn_m   ! Mean d(qns)/dT over all categories
217      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_devap_m ! Mean d(evap)/dT over all categories
218      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zalb_m    ! Mean albedo over all categories
219      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztem_m    ! Mean temperature over all categories
220      !!----------------------------------------------------------------------
221      !
222      WHERE ( at_i (:,:) > 0._wp )   ; z1_at_i(:,:) = 1._wp / at_i (:,:)
223      ELSEWHERE                      ; z1_at_i(:,:) = 0._wp
224      END WHERE
225     
226      SELECT CASE( k_flxdist )       !==  averaged on all ice categories  ==!
227      !
228      CASE( 0 , 1 )
229         !
230         ALLOCATE( z_qns_m(jpi,jpj), z_qsr_m(jpi,jpj), z_dqn_m(jpi,jpj), z_evap_m(jpi,jpj), z_devap_m(jpi,jpj) ) 
231         !
232         z_qns_m  (:,:) = SUM( a_i(:,:,:) * pqns_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
233         z_qsr_m  (:,:) = SUM( a_i(:,:,:) * pqsr_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
234         z_dqn_m  (:,:) = SUM( a_i(:,:,:) * pdqn_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
235         z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
236         z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
237         DO jl = 1, jpl
238            pqns_ice  (:,:,jl) = z_qns_m (:,:)
239            pqsr_ice  (:,:,jl) = z_qsr_m (:,:)
240            pdqn_ice  (:,:,jl) = z_dqn_m  (:,:)
241            pevap_ice (:,:,jl) = z_evap_m(:,:)
242            pdevap_ice(:,:,jl) = z_devap_m(:,:)
243         END DO
244         !
245         DEALLOCATE( z_qns_m, z_qsr_m, z_dqn_m, z_evap_m, z_devap_m ) 
246         !
247      END SELECT
248      !
249      SELECT CASE( k_flxdist )       !==  redistribution on all ice categories  ==!
250      !
251      CASE( 1 , 2 )
252         !
253         ALLOCATE( zalb_m(jpi,jpj), ztem_m(jpi,jpj) ) 
254         !
255         zalb_m(:,:) = SUM( a_i(:,:,:) * palb_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
256         ztem_m(:,:) = SUM( a_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
257         DO jl = 1, jpl
258            pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice  (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
259            pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
260            pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
261         END DO
262         !
263         DEALLOCATE( zalb_m, ztem_m ) 
264         !
265      END SELECT
266      !
267   END SUBROUTINE ice_flx_dist
268
269
270   SUBROUTINE ice_sbc_init
271      !!-------------------------------------------------------------------
272      !!                  ***  ROUTINE ice_sbc_init  ***
273      !!
274      !! ** Purpose :   Physical constants and parameters linked to the ice dynamics
275      !!     
276      !! ** Method  :   Read the namsbc namelist and check the ice-dynamic
277      !!              parameter values called at the first timestep (nit000)
278      !!
279      !! ** input   :   Namelist namsbc
280      !!-------------------------------------------------------------------
281      INTEGER ::   ios, ioptio   ! Local integer
282      !!
283      NAMELIST/namsbc/ rn_cio, rn_blow_s, nn_flxdist, ln_cndflx, ln_cndemulate
284      !!-------------------------------------------------------------------
285      !
286      READ  ( numnam_ice_ref, namsbc, IOSTAT = ios, ERR = 901)
287901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc in reference namelist' )
288      READ  ( numnam_ice_cfg, namsbc, IOSTAT = ios, ERR = 902 )
289902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc in configuration namelist' )
290      IF(lwm) WRITE( numoni, namsbc )
291      !
292      IF(lwp) THEN                     ! control print
293         WRITE(numout,*)
294         WRITE(numout,*) 'ice_sbc_init: ice parameters for ice dynamics '
295         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
296         WRITE(numout,*) '   Namelist namsbc:'
297         WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio        = ', rn_cio
298         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s     = ', rn_blow_s
299         WRITE(numout,*) '      Multicategory heat flux formulation              nn_flxdist    = ', nn_flxdist
300         WRITE(numout,*) '      Use conduction flux as surface condition         ln_cndflx     = ', ln_cndflx
301         WRITE(numout,*) '         emulate conduction flux                       ln_cndemulate = ', ln_cndemulate
302      ENDIF
303      !
304      IF(lwp) WRITE(numout,*)
305      SELECT CASE( nn_flxdist )         ! SI3 Multi-category heat flux formulation
306      CASE( -1  )
307         IF(lwp) WRITE(numout,*) '   SI3: use per-category fluxes (nn_flxdist = -1) '
308      CASE(  0  )
309         IF(lwp) WRITE(numout,*) '   SI3: use average per-category fluxes (nn_flxdist = 0) '
310      CASE(  1  )
311         IF(lwp) WRITE(numout,*) '   SI3: use average then redistribute per-category fluxes (nn_flxdist = 1) '
312         IF( ln_cpl )         CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for SI3 in coupled mode must be /=1' )
313      CASE(  2  )
314         IF(lwp) WRITE(numout,*) '   SI3: Redistribute a single flux over categories (nn_flxdist = 2) '
315         IF( .NOT. ln_cpl )   CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for SI3 in forced mode must be /=2' )
316      CASE DEFAULT
317         CALL ctl_stop( 'ice_thd_init: SI3 option, nn_flxdist, should be between -1 and 2' )
318      END SELECT
319      !
320   END SUBROUTINE ice_sbc_init
321
322#else
323   !!----------------------------------------------------------------------
324   !!   Default option :         Empty module         NO SI3 sea-ice model
325   !!----------------------------------------------------------------------
326#endif
327
328   !!======================================================================
329END MODULE icesbc
Note: See TracBrowser for help on using the repository browser.