source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icesbc.F90 @ 13140

Last change on this file since 13140 was 12546, checked in by orioltp, 12 months ago

Adding precision specification in hardcoded reals and other modifications to allow compilation without forcing reals without precision specification to a certain value through compiler flags

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