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.
iceforcing.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

Last change on this file since 8414 was 8414, checked in by clem, 7 years ago

changing names
File size: 13.9 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 icealbedo       ! 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      SELECT CASE( ksbc )
67         CASE( jp_usr     )   ;    CALL usrdef_sbc_ice_tau( kt )                 ! user defined formulation
68         CASE( jp_blk     )   ;    CALL blk_ice_tau                              ! Bulk formulation
69         CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled   formulation
70      END SELECT
71
72      IF( ln_mixcpl) THEN                                                        ! Case of a mixed Bulk/Coupled formulation
73                                   CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
74         DO jj = 2, jpjm1
75            DO ji = 2, jpim1
76               utau_ice(ji,jj) = utau_ice(ji,jj) * xcplmask(ji,jj,0) + zutau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
77               vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
78            END DO
79         END DO
80         CALL lbc_lnk( utau_ice, 'U', -1. )
81         CALL lbc_lnk( vtau_ice, 'V', -1. )
82      ENDIF
83
84      IF( nn_timing == 1 )   CALL timing_stop('ice_forcing_tau')
85      !
86   END SUBROUTINE ice_forcing_tau
87
88   
89   SUBROUTINE ice_forcing_flx( kt, ksbc )
90      !!---------------------------------------------------------------------
91      !!                  ***  ROUTINE ice_forcing_flx  ***
92      !!
93      !! ** Purpose : provide surface boundary condition for sea ice (flux)
94      !!
95      !! ** Action  : It provides the following fields used in sea ice model:
96      !!                fr1_i0  , fr2_i0                         = 1sr & 2nd fraction of qsr penetration in ice  [%]
97      !!                emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s]
98      !!                sprecip                                  = solid precipitation                           [Kg/m2/s]
99      !!                evap_ice                                 = sublimation                                   [Kg/m2/s]
100      !!                qsr_tot , qns_tot                        = solar & non solar heat flux (total)           [W/m2]
101      !!                qsr_ice , qns_ice                        = solar & non solar heat flux over ice          [W/m2]
102      !!                dqns_ice                                 = non solar  heat sensistivity                  [W/m2]
103      !!                qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2]
104      !!            + some fields that are not used outside this module:
105      !!                qla_ice                                  = latent heat flux over ice                     [W/m2]
106      !!                dqla_ice                                 = latent heat sensistivity                      [W/m2]
107      !!                tprecip                                  = total  precipitation                          [Kg/m2/s]
108      !!                alb_ice                                  = albedo above sea ice
109      !!---------------------------------------------------------------------
110      INTEGER, INTENT(in) ::   kt      ! ocean time step
111      INTEGER, INTENT(in) ::   ksbc    ! type of sbc flux ( 1 = user defined formulation,
112                                       !                    3 = bulk formulation,
113                                       !                    4 = Pure Coupled formulation)
114      INTEGER  ::   ji, jj, jl                                ! dummy loop index
115      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
116      REAL(wp), DIMENSION(jpi,jpj)     ::   zalb              ! 2D workspace
117      !!----------------------------------------------------------------------
118
119      IF( nn_timing == 1 )   CALL timing_start('ice_forcing_flx')
120
121      ! --- cloud-sky and overcast-sky ice albedos --- !
122      CALL ice_albedo( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os )
123
124      ! albedo depends on cloud fraction because of non-linear spectral effects
125      DO jl = 1, jpl
126         DO jj = 2, jpjm1
127            DO ji = 2, jpim1                               
128               alb_ice(ji,jj,jl) = ( 1. - cldf_ice ) * zalb_cs(ji,jj,jl) + cldf_ice * zalb_os(ji,jj,jl)
129            END DO
130         END DO
131      END DO
132      CALL lbc_lnk( alb_ice, 'T', 1. )
133     
134      ! --- fluxes over sea ice --- !
135      SELECT CASE( ksbc )
136         CASE( jp_usr )                                          ! user defined formulation
137                                   CALL usrdef_sbc_ice_flx( kt )
138
139         CASE( jp_blk )                                          ! bulk formulation
140                                   CALL blk_ice_flx( t_su, alb_ice )
141            IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
142            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
143
144         CASE ( jp_purecpl )                                     ! coupled formulation
145                                   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
146            IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
147      END SELECT
148
149      ! --- albedo output --- !
150      WHERE( at_i_b <= epsi06 )  ;  zalb(:,:) = rn_alb_oce
151      ELSEWHERE                  ;  zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
152      END WHERE
153      IF( iom_use('icealb'  ) )  CALL iom_put( "icealb" , zalb(:,:) )   ! ice albedo
154
155      zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b )
156      IF( iom_use('albedo'  ) )  CALL iom_put( "albedo" , zalb(:,:) )   ! surface albedo
157      !
158      !
159      IF( nn_timing == 1 )   CALL timing_stop('ice_forcing_flx')
160
161   END SUBROUTINE ice_forcing_flx
162
163
164   SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
165      !!---------------------------------------------------------------------
166      !!                  ***  ROUTINE ice_lim_flx  ***
167      !!
168      !! ** Purpose :   update the ice surface boundary condition by averaging and / or
169      !!                redistributing fluxes on ice categories
170      !!
171      !! ** Method  :   average then redistribute
172      !!
173      !! ** Action  :
174      !!---------------------------------------------------------------------
175      INTEGER                   , INTENT(in   ) ::   k_limflx   ! =-1 do nothing; =0 average ;
176      !                                                         ! = 1 average and redistribute ; =2 redistribute
177      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature
178      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo
179      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqns_ice   ! non solar flux
180      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux
181      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdqn_ice   ! non solar flux sensitivity
182      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pevap_ice  ! sublimation
183      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity
184      !
185      INTEGER  ::   jl      ! dummy loop index
186      !
187      REAL(wp), DIMENSION(jpi,jpj) :: zalb_m    ! Mean albedo over all categories
188      REAL(wp), DIMENSION(jpi,jpj) :: ztem_m    ! Mean temperature over all categories
189      !
190      REAL(wp), DIMENSION(jpi,jpj) :: z_qsr_m   ! Mean solar heat flux over all categories
191      REAL(wp), DIMENSION(jpi,jpj) :: z_qns_m   ! Mean non solar heat flux over all categories
192      REAL(wp), DIMENSION(jpi,jpj) :: z_evap_m  ! Mean sublimation over all categories
193      REAL(wp), DIMENSION(jpi,jpj) :: z_dqn_m   ! Mean d(qns)/dT over all categories
194      REAL(wp), DIMENSION(jpi,jpj) :: z_devap_m ! Mean d(evap)/dT over all categories
195      !!----------------------------------------------------------------------
196      !
197      IF( nn_timing == 1 )  CALL timing_start('ice_lim_flx')
198      !
199      SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==!
200      CASE( 0 , 1 )
201         !
202         z_qns_m  (:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
203         z_qsr_m  (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
204         z_dqn_m  (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
205         z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) )
206         z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) )
207         DO jl = 1, jpl
208            pdqn_ice  (:,:,jl) = z_dqn_m(:,:)
209            pdevap_ice(:,:,jl) = z_devap_m(:,:)
210         END DO
211         !
212         DO jl = 1, jpl
213            pqns_ice (:,:,jl) = z_qns_m(:,:)
214            pqsr_ice (:,:,jl) = z_qsr_m(:,:)
215            pevap_ice(:,:,jl) = z_evap_m(:,:)
216         END DO
217         !
218      END SELECT
219      !
220      SELECT CASE( k_limflx )                              !==  redistribution on all ice categories  ==!
221      CASE( 1 , 2 )
222         !
223         zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) )
224         ztem_m(:,:) = fice_ice_ave ( ptn_ice  (:,:,:) )
225         DO jl = 1, jpl
226            pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice  (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
227            pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
228            pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
229         END DO
230         !
231      END SELECT
232      !
233      IF( nn_timing == 1 )  CALL timing_stop('ice_lim_flx')
234      !
235   END SUBROUTINE ice_lim_flx
236
237
238   FUNCTION fice_cell_ave ( ptab )
239      !!--------------------------------------------------------------------------
240      !! * Compute average over categories, for grid cell (ice covered and free ocean)
241      !!--------------------------------------------------------------------------
242      REAL(wp), DIMENSION(jpi,jpj) :: fice_cell_ave
243      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT (in) :: ptab
244      INTEGER :: jl
245      !!--------------------------------------------------------------------------
246
247      fice_cell_ave (:,:) = 0._wp
248      DO jl = 1, jpl
249         fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl)
250      END DO
251
252   END FUNCTION fice_cell_ave
253
254
255   FUNCTION fice_ice_ave ( ptab )
256      !!--------------------------------------------------------------------------
257      !! * Compute average over categories, for ice covered part of grid cell
258      !!--------------------------------------------------------------------------
259      REAL(wp), DIMENSION(jpi,jpj) :: fice_ice_ave
260      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in) :: ptab
261     !!--------------------------------------------------------------------------
262
263      WHERE ( at_i (:,:) > 0.0_wp ) ; fice_ice_ave (:,:) = fice_cell_ave( ptab (:,:,:) ) / at_i (:,:)
264      ELSEWHERE                     ; fice_ice_ave (:,:) = 0.0_wp
265      END WHERE
266     
267   END FUNCTION fice_ice_ave
268
269#endif
270
271   !!======================================================================
272END MODULE iceforcing
Note: See TracBrowser for help on using the repository browser.