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.
usrdef_sbc.F90 in NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/tests/ICE_RHEO/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90 @ 14063

Last change on this file since 14063 was 14063, checked in by laurent, 3 years ago

Catch up with trunk at r14060
File size: 10.9 KB
Line 
1MODULE usrdef_sbc
2   !!======================================================================
3   !!                       ***  MODULE  usrdef_sbc  ***
4   !!
5   !!                      ===  ICE_RHEO configuration  ===
6   !!
7   !! User defined :   surface forcing of a user configuration
8   !!======================================================================
9   !! History :  4.0   ! 2016-03  (S. Flavoni, G. Madec)  user defined interface
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   usr_def_sbc    : user defined surface bounday conditions in ICE_RHEO case
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE sbc_oce         ! Surface boundary condition: ocean fields
18   USE sbc_ice         ! Surface boundary condition: ice fields
19   USE phycst          ! physical constants
20   USE ice, ONLY       : at_i_b, a_i_b, at_i, u_ice, v_ice
21   USE icethd_dh       ! for CALL ice_thd_snwblow
22   USE sbc_phy, ONLY : pp_cldf
23   !
24   USE in_out_manager  ! I/O manager
25   USE lib_mpp         ! distribued memory computing library
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   usrdef_sbc_oce      ! routine called by sbcmod.F90 for sbc ocean
33   PUBLIC   usrdef_sbc_ice_tau  ! routine called by icestp.F90 for ice dynamics
34   PUBLIC   usrdef_sbc_ice_flx  ! routine called by icestp.F90 for ice thermo
35
36   !! * Substitutions
37#  include "do_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
40   !! $Id: usrdef_sbc.F90 10074 2018-08-28 16:15:49Z nicolasmartin $
41   !! Software governed by the CeCILL license (see ./LICENSE)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE usrdef_sbc_oce( kt, Kbb )
46      !!---------------------------------------------------------------------
47      !!                    ***  ROUTINE usr_def_sbc  ***
48      !!             
49      !! ** Purpose :   provide at each time-step the surface boundary
50      !!              condition, i.e. the momentum, heat and freshwater fluxes.
51      !!
52      !! ** Method  :   all 0 fields, for ICE_RHEO case
53      !!                CAUTION : never mask the surface stress field !
54      !!
55      !! ** Action  : - set to ZERO all the ocean surface boundary condition, i.e.   
56      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
57      !!
58      !!----------------------------------------------------------------------
59      INTEGER, INTENT(in) ::   kt   ! ocean time step
60      INTEGER, INTENT(in) ::   Kbb  ! ocean time index
61      INTEGER ::    ij0, ij1, ii0, ii1, jj, ji   ! loop indices
62      REAL(wp) ::   zrhoco          ! ocean density and drag coefficient product
63      !!---------------------------------------------------------------------
64      !
65      IF( kt == nit000 ) THEN
66         !
67         !IF(lwp)   WRITE(numout,*)' usrdef_sbc_oce : ICE_RHEO case: ocean boudary conditions'
68
69         utau(:,:) = 0._wp
70         utau(:,:) = 0._wp
71
72         !ij0 = 1   ;   ij1 = 25                       ! set boundary condition
73         !ii0 = 975   ;   ii1 = 1000
74         !DO jj = mj0(ij0), mj1(ij1)
75         !   DO ji = mi0(ii0), mi1(ii1)
76         !      utau(ji,jj) = -utau_ice(ji,jj)
77         !      vtau(ji,jj) = -vtau_ice(ji,jj)
78         !   END DO
79         !END DO
80
81         taum(:,:) = 0._wp   ! assume these are not used
82         wndm(:,:) = 0._wp
83         !
84         emp (:,:) = 0._wp
85         sfx (:,:) = 0._wp
86         qns (:,:) = 0._wp
87         qsr (:,:) = 0._wp
88         !
89         utau_b(:,:) = 0._wp 
90         vtau_b(:,:) = 0._wp
91         emp_b (:,:) = 0._wp
92         sfx_b (:,:) = 0._wp
93         qns_b (:,:) = 0._wp
94         !
95      ENDIF
96      !
97   END SUBROUTINE usrdef_sbc_oce
98
99   SUBROUTINE usrdef_sbc_ice_tau( kt )
100      !!---------------------------------------------------------------------
101      !!                     ***  ROUTINE usrdef_sbc_ice_tau  ***
102      !!
103      !! ** Purpose :   provide the surface boundary (momentum) condition over
104      !sea-ice
105      !!---------------------------------------------------------------------
106      INTEGER, INTENT(in) ::   kt   ! ocean time step
107      INTEGER ::   jj, ji   ! loop indices
108
109      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f       ! relative wind module and components at F-point
110      REAL(wp) ::   zwndi_t , zwndj_t                 ! relative wind components at T-point
111      REAL(wp), DIMENSION(jpi,jpj) ::   windu, windv  ! wind components (idealised forcing)
112      REAL(wp), PARAMETER ::   r_vfac = 1._wp         ! relative velocity (make 0 for absolute velocity)
113      REAL(wp), PARAMETER ::   Rwind = -0.8_wp        ! ratio of wind components
114      REAL(wp), PARAMETER ::   Umax = 15._wp          ! maximum wind speed (m/s)
115      REAL(wp), PARAMETER ::   d = 2000._wp           ! size of the domain (km)
116      REAL(wp), PARAMETER ::   res = 2._wp            ! gridcell size
117      REAL(wp), PARAMETER ::   zrhoa  = 1.22          ! Air density kg/m3
118      REAL(wp), PARAMETER ::   Cd_atm =  1.4e-3       ! transfer coefficient over ice
119      !!---------------------------------------------------------------------
120      ! extra code for test case
121      IF( kt==nit000 .AND. lwp)   WRITE(numout,*)' usrdef_sbc_ice : ICE_RHEO case: analytical stress forcing'
122
123      DO_2D( 0, 0, 0, 0 )
124         ! wind spins up over 6 hours, factor 1000 to balance the units
125         windu(ji,jj) = Umax/sqrt(d*1000)*(d-2*mig(ji)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*min(kt*30./21600,1.)
126         windv(ji,jj) = Umax/sqrt(d*1000)*(d-2*mjg(jj)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*Rwind*min(kt*30./21600,1.)
127      END_2D
128      CALL lbc_lnk_multi( 'usrdef_sbc', windu, 'U', -1., windv, 'V', -1. )
129
130      wndm_ice(:,:) = 0._wp      !!gm brutal....
131
132      ! ------------------------------------------------------------ !
133      !    Wind module relative to the moving ice ( U10m - U_ice )   !
134      ! ------------------------------------------------------------ !
135      ! C-grid ice dynamics :   U & V-points (same as ocean)
136      DO_2D( 0, 0, 0, 0 )
137         zwndi_t = (  windu(ji,jj) - r_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
138         zwndj_t = (  windv(ji,jj) - r_vfac * 0.5 * ( v_ice(ji,jj-1) + v_ice(ji,jj) )  )
139         wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
140      END_2D
141      CALL lbc_lnk( 'usrdef_sbc', wndm_ice, 'T',  1. )
142
143      !!gm brutal....
144      utau_ice  (:,:) = 0._wp
145      vtau_ice  (:,:) = 0._wp
146      !!gm end
147
148      ! ------------------------------------------------------------ !
149      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
150      ! ------------------------------------------------------------ !
151      ! C-grid ice dynamics :   U & V-points (same as ocean)
152      DO_2D( 0, 0, 0, 0 )
153         utau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
154            &          * ( 0.5 * (windu(ji+1,jj) + windu(ji,jj) ) - r_vfac * u_ice(ji,jj) )
155         vtau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
156            &          * ( 0.5 * (windv(ji,jj+1) + windv(ji,jj) ) - r_vfac * v_ice(ji,jj) )
157      END_2D
158      CALL lbc_lnk_multi( 'usrdef_sbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
159      !
160   END SUBROUTINE usrdef_sbc_ice_tau
161
162   SUBROUTINE usrdef_sbc_ice_flx( kt, phs, phi )
163      !!---------------------------------------------------------------------
164      !!                     ***  ROUTINE usrdef_sbc_ice_flx  ***
165      !!
166      !! ** Purpose :   provide the surface boundary (flux) condition over sea-ice
167      !!---------------------------------------------------------------------
168      INTEGER, INTENT(in) ::   kt   ! ocean time step
169      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phs    ! snow thickness
170      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi    ! ice thickness
171      !!
172      REAL(wp) ::   zfr1, zfr2                 ! local variables
173      REAL(wp), DIMENSION(jpi,jpj) ::   zsnw   ! snw distribution after wind blowing
174      !!---------------------------------------------------------------------
175      !
176      IF( kt==nit000 .AND. lwp)   WRITE(numout,*)' usrdef_sbc_ice : ICE_RHEO case: NO flux forcing'
177      !
178      ! ocean variables (renaming)
179      emp_oce (:,:)   = 0._wp   ! uniform value for freshwater budget (E-P)
180      qsr_oce (:,:)   = 0._wp   ! uniform value for     solar radiation
181      qns_oce (:,:)   = 0._wp   ! uniform value for non-solar radiation
182
183      ! ice variables
184      alb_ice (:,:,:) = 0.7_wp  ! useless
185      qsr_ice (:,:,:) = 0._wp   ! uniform value for     solar radiation
186      qns_ice (:,:,:) = 0._wp   ! uniform value for non-solar radiation
187      sprecip (:,:)   = 0._wp   ! uniform value for snow precip
188      evap_ice(:,:,:) = 0._wp   ! uniform value for sublimation
189
190      ! ice fields deduced from above
191      zsnw(:,:) = 1._wp
192      !!CALL lim_thd_snwblow( at_i_b, zsnw )  ! snow distribution over ice after wind blowing
193      emp_ice  (:,:)   = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw(:,:)
194      emp_oce  (:,:)   = emp_oce(:,:) - sprecip(:,:) * (1._wp - zsnw(:,:) )
195      qevap_ice(:,:,:) =   0._wp
196      qprec_ice(:,:)   =   rhos * ( sst_m(:,:) * rcpi - rLfus ) * tmask(:,:,1) !  in J/m3
197      qemp_oce (:,:)   = - emp_oce(:,:) * sst_m(:,:) * rcp
198      qemp_ice (:,:)   =   sprecip(:,:) * zsnw * ( sst_m(:,:) * rcpi - rLfus ) * tmask(:,:,1) ! solid precip (only)
199
200      ! total fluxes
201      emp_tot (:,:) = emp_ice  + emp_oce
202      qns_tot (:,:) = at_i_b(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:)
203      qsr_tot (:,:) = at_i_b(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
204
205      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
206      zfr1 = ( 0.18 * ( 1.0 - pp_cldf ) + 0.35 * pp_cldf )            ! transmission when hi>10cm
207      zfr2 = ( 0.82 * ( 1.0 - pp_cldf ) + 0.65 * pp_cldf )            ! zfr2 such that zfr1 + zfr2 to equal 1
208      !
209      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm 
210         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
211      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
212         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
213      ELSEWHERE                                                         ! zero when hs>0
214         qtr_ice_top(:,:,:) = 0._wp 
215      END WHERE
216         
217   END SUBROUTINE usrdef_sbc_ice_flx
218
219
220   !!======================================================================
221END MODULE usrdef_sbc
Note: See TracBrowser for help on using the repository browser.