source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90 @ 8426

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

last routine names to be changed

File size: 14.2 KB
Line 
1MODULE iceforcing
2   !!======================================================================
3   !!                       ***  MODULE  iceforcing  ***
4   !! call surface forcing for sea ice model
5   !!=====================================================================
6   !! History :  4.0  ! 2017-08  (C. Rousset) Original code
7   !!----------------------------------------------------------------------
8#if defined key_lim3
9   !!----------------------------------------------------------------------
10   !!   'key_lim3' :                                  LIM 3.0 sea-ice model
11   !!----------------------------------------------------------------------
12   USE oce             ! ocean dynamics and tracers
13   USE dom_oce         ! ocean space and time domain
14   USE ice             ! LIM-3: ice variables
15   !
16   USE sbc_oce         ! Surface boundary condition: ocean fields
17   USE sbc_ice         ! Surface boundary condition: ice   fields
18   USE usrdef_sbc      ! user defined: surface boundary condition
19   USE sbcblk          ! Surface boundary condition: bulk
20   USE sbccpl          ! Surface boundary condition: coupled interface
21   USE icealb          ! ice albedo
22   !
23   USE iom             ! I/O manager library
24   USE in_out_manager  ! I/O manager
25   USE lbclnk          ! lateral boundary condition - MPP link
26   USE lib_mpp         ! MPP library
27   USE lib_fortran     !
28   USE timing          ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC ice_forcing_tau  ! routine called by icestp.F90
34   PUBLIC ice_forcing_flx  ! routine called by icestp.F90
35
36   !! * Substitutions
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 4.0 , UCL NEMO Consortium (2011)
40   !! $Id: icestp.F90 8319 2017-07-11 15:00:44Z clem $
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice )
46      !!---------------------------------------------------------------------
47      !!                  ***  ROUTINE ice_forcing_tau  ***
48      !!
49      !! ** Purpose : provide surface boundary condition for sea ice (momentum)
50      !!
51      !! ** Action  : It provides the following fields:
52      !!              utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2]
53      !!---------------------------------------------------------------------
54      INTEGER, INTENT(in) ::   kt      ! ocean time step
55      INTEGER, INTENT(in) ::   ksbc    ! type of sbc flux ( 1 = user defined formulation,
56                                       !                    3 = bulk formulation,
57                                       !                    4 = Pure Coupled formulation)
58      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) ::   utau_ice, vtau_ice 
59      !!
60      INTEGER  ::   ji, jj                 ! dummy loop index
61      REAL(wp), DIMENSION(jpi,jpj) ::   zutau_ice, zvtau_ice 
62      !!----------------------------------------------------------------------
63
64      IF( nn_timing == 1 )   CALL timing_start('ice_forcing_tau')
65
66      IF( kt == nit000 .AND. lwp ) THEN
67         WRITE(numout,*)
68         WRITE(numout,*)'ice_forcing_tau'
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_tau                              ! Bulk formulation
75         CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled   formulation
76      END SELECT
77
78      IF( ln_mixcpl) THEN                                                        ! Case of a mixed Bulk/Coupled formulation
79                                   CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
80         DO jj = 2, jpjm1
81            DO ji = 2, jpim1
82               utau_ice(ji,jj) = utau_ice(ji,jj) * xcplmask(ji,jj,0) + zutau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
83               vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
84            END DO
85         END DO
86         CALL lbc_lnk_multi( utau_ice, 'U', -1., vtau_ice, 'V', -1. )
87      ENDIF
88
89      IF( nn_timing == 1 )   CALL timing_stop('ice_forcing_tau')
90      !
91   END SUBROUTINE ice_forcing_tau
92
93   
94   SUBROUTINE ice_forcing_flx( kt, ksbc )
95      !!---------------------------------------------------------------------
96      !!                  ***  ROUTINE ice_forcing_flx  ***
97      !!
98      !! ** Purpose : provide surface boundary condition for sea ice (flux)
99      !!
100      !! ** Action  : It provides the following fields used in sea ice model:
101      !!                fr1_i0  , fr2_i0                         = 1sr & 2nd fraction of qsr penetration in ice  [%]
102      !!                emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s]
103      !!                sprecip                                  = solid precipitation                           [Kg/m2/s]
104      !!                evap_ice                                 = sublimation                                   [Kg/m2/s]
105      !!                qsr_tot , qns_tot                        = solar & non solar heat flux (total)           [W/m2]
106      !!                qsr_ice , qns_ice                        = solar & non solar heat flux over ice          [W/m2]
107      !!                dqns_ice                                 = non solar  heat sensistivity                  [W/m2]
108      !!                qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2]
109      !!            + some fields that are not used outside this module:
110      !!                qla_ice                                  = latent heat flux over ice                     [W/m2]
111      !!                dqla_ice                                 = latent heat sensistivity                      [W/m2]
112      !!                tprecip                                  = total  precipitation                          [Kg/m2/s]
113      !!                alb_ice                                  = albedo above sea ice
114      !!---------------------------------------------------------------------
115      INTEGER, INTENT(in) ::   kt      ! ocean time step
116      INTEGER, INTENT(in) ::   ksbc    ! type of sbc flux ( 1 = user defined formulation,
117                                       !                    3 = bulk formulation,
118                                       !                    4 = Pure Coupled formulation)
119      INTEGER  ::   ji, jj, jl                                ! dummy loop index
120      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
121      REAL(wp), DIMENSION(jpi,jpj)     ::   zalb              ! 2D workspace
122      !!----------------------------------------------------------------------
123
124      IF( nn_timing == 1 )   CALL timing_start('ice_forcing_flx')
125
126      IF( kt == nit000 .AND. lwp ) THEN
127         WRITE(numout,*)
128         WRITE(numout,*)'ice_forcing_flx'
129         WRITE(numout,*)'~~~~~~~~~~~~~~~'
130      ENDIF
131
132      ! --- cloud-sky and overcast-sky ice albedos --- !
133      CALL ice_alb( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os )
134
135      ! albedo depends on cloud fraction because of non-linear spectral effects
136      DO jl = 1, jpl
137         DO jj = 2, jpjm1
138            DO ji = 2, jpim1                               
139               alb_ice(ji,jj,jl) = ( 1. - cldf_ice ) * zalb_cs(ji,jj,jl) + cldf_ice * zalb_os(ji,jj,jl)
140            END DO
141         END DO
142      END DO
143      CALL lbc_lnk( alb_ice, 'T', 1. )
144     
145      ! --- fluxes over sea ice --- !
146      SELECT CASE( ksbc )
147         CASE( jp_usr )                                          ! user defined formulation
148                                   CALL usrdef_sbc_ice_flx( kt )
149
150         CASE( jp_blk )                                          ! bulk formulation
151                                   CALL blk_ice_flx( t_su, alb_ice )
152            IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
153            IF( nn_limflx /= 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
154
155         CASE ( jp_purecpl )                                     ! coupled formulation
156                                   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
157            IF( nn_limflx == 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
158      END SELECT
159
160      ! --- albedo output --- !
161      WHERE( at_i_b <= epsi06 )  ;  zalb(:,:) = rn_alb_oce
162      ELSEWHERE                  ;  zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
163      END WHERE
164      IF( iom_use('icealb'  ) )  CALL iom_put( "icealb" , zalb(:,:) )   ! ice albedo
165
166      zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b )
167      IF( iom_use('albedo'  ) )  CALL iom_put( "albedo" , zalb(:,:) )   ! surface albedo
168      !
169      !
170      IF( nn_timing == 1 )   CALL timing_stop('ice_forcing_flx')
171
172   END SUBROUTINE ice_forcing_flx
173
174
175   SUBROUTINE ice_flx_dist( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
176      !!---------------------------------------------------------------------
177      !!                  ***  ROUTINE ice_flx_dist  ***
178      !!
179      !! ** Purpose :   update the ice surface boundary condition by averaging and / or
180      !!                redistributing fluxes on ice categories
181      !!
182      !! ** Method  :   average then redistribute
183      !!
184      !! ** Action  :
185      !!---------------------------------------------------------------------
186      INTEGER                   , INTENT(in   ) ::   k_limflx   ! =-1 do nothing; =0 average ;
187      !                                                         ! = 1 average and redistribute ; =2 redistribute
188      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature
189      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo
190      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqns_ice   ! non solar flux
191      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux
192      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdqn_ice   ! non solar flux sensitivity
193      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pevap_ice  ! sublimation
194      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity
195      !
196      INTEGER  ::   jl      ! dummy loop index
197      !
198      REAL(wp), DIMENSION(jpi,jpj) :: zalb_m    ! Mean albedo over all categories
199      REAL(wp), DIMENSION(jpi,jpj) :: ztem_m    ! Mean temperature over all categories
200      !
201      REAL(wp), DIMENSION(jpi,jpj) :: z_qsr_m   ! Mean solar heat flux over all categories
202      REAL(wp), DIMENSION(jpi,jpj) :: z_qns_m   ! Mean non solar heat flux over all categories
203      REAL(wp), DIMENSION(jpi,jpj) :: z_evap_m  ! Mean sublimation over all categories
204      REAL(wp), DIMENSION(jpi,jpj) :: z_dqn_m   ! Mean d(qns)/dT over all categories
205      REAL(wp), DIMENSION(jpi,jpj) :: z_devap_m ! Mean d(evap)/dT over all categories
206      !!----------------------------------------------------------------------
207      !
208      IF( nn_timing == 1 )  CALL timing_start('ice_flx_dist')
209      !
210      SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==!
211      CASE( 0 , 1 )
212         !
213         z_qns_m  (:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
214         z_qsr_m  (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
215         z_dqn_m  (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
216         z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) )
217         z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) )
218         DO jl = 1, jpl
219            pdqn_ice  (:,:,jl) = z_dqn_m(:,:)
220            pdevap_ice(:,:,jl) = z_devap_m(:,:)
221         END DO
222         !
223         DO jl = 1, jpl
224            pqns_ice (:,:,jl) = z_qns_m(:,:)
225            pqsr_ice (:,:,jl) = z_qsr_m(:,:)
226            pevap_ice(:,:,jl) = z_evap_m(:,:)
227         END DO
228         !
229      END SELECT
230      !
231      SELECT CASE( k_limflx )                              !==  redistribution on all ice categories  ==!
232      CASE( 1 , 2 )
233         !
234         zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) )
235         ztem_m(:,:) = fice_ice_ave ( ptn_ice  (:,:,:) )
236         DO jl = 1, jpl
237            pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice  (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
238            pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
239            pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
240         END DO
241         !
242      END SELECT
243      !
244      IF( nn_timing == 1 )  CALL timing_stop('ice_flx_dist')
245      !
246   END SUBROUTINE ice_flx_dist
247
248
249   FUNCTION fice_cell_ave ( ptab )
250      !!--------------------------------------------------------------------------
251      !! * Compute average over categories, for grid cell (ice covered and free ocean)
252      !!--------------------------------------------------------------------------
253      REAL(wp), DIMENSION(jpi,jpj) :: fice_cell_ave
254      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT (in) :: ptab
255      INTEGER :: jl
256      !!--------------------------------------------------------------------------
257
258      fice_cell_ave (:,:) = 0._wp
259      DO jl = 1, jpl
260         fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl)
261      END DO
262
263   END FUNCTION fice_cell_ave
264
265
266   FUNCTION fice_ice_ave ( ptab )
267      !!--------------------------------------------------------------------------
268      !! * Compute average over categories, for ice covered part of grid cell
269      !!--------------------------------------------------------------------------
270      REAL(wp), DIMENSION(jpi,jpj) :: fice_ice_ave
271      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in) :: ptab
272     !!--------------------------------------------------------------------------
273
274      WHERE ( at_i (:,:) > 0.0_wp ) ; fice_ice_ave (:,:) = fice_cell_ave( ptab (:,:,:) ) / at_i (:,:)
275      ELSEWHERE                     ; fice_ice_ave (:,:) = 0.0_wp
276      END WHERE
277     
278   END FUNCTION fice_ice_ave
279
280#endif
281
282   !!======================================================================
283END MODULE iceforcing
Note: See TracBrowser for help on using the repository browser.