source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90 @ 11322

Last change on this file since 11322 was 11322, checked in by flemarie, 2 years ago

First implementation of ABL (see ticket #2131)

  • Update reference and cfg namelists for ORCA2_ICE_ABL
  • Run ABL over the ocean and BLK over sea-ice (ABL over sea-ice to come)
  • Bug fix in computation of pblh (+ add option to smooth pblh)
File size: 18.0 KB
Line 
1MODULE sbcabl
2   !!======================================================================
3   !!                       ***  MODULE  sbcabl  ***
4   !! Ocean forcing:  momentum, heat and freshwater flux formulation
5   !!                 derived from an ABL model
6   !!=====================================================================
7   !! History :  4.0  !  2019-03  (F. Lemarié & G. Samson)  Original code
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   sbc_abl_init  : Initialization of ABL model based on namelist options
12   !!   sbc_abl       : driver for the computation of momentum, heat and freshwater
13   !!                   fluxes over ocean via the ABL model
14   !!----------------------------------------------------------------------
15   USE abl            ! ABL
16   USE par_abl        ! abl parameters
17   USE ablmod
18
19   USE phycst         ! physical constants
20   USE fldread        ! read input fields
21   USE sbc_oce        ! Surface boundary condition: ocean fields
22   USE sbcblk         ! Surface boundary condition: bulk formulae
23   USE dom_oce, ONLY  : tmask
24   !
25   USE iom            ! I/O manager library
26   USE in_out_manager ! I/O manager
27   USE lib_mpp        ! distribued memory computing library
28   USE lib_fortran    ! to use key_nosignedzero
29   USE timing         ! Timing
30   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
31   USE prtctl         ! Print control
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   sbc_abl_init       ! routine called in sbcmod module
37   PUBLIC   sbc_abl            ! routine called in sbcmod module
38
39   !! * Substitutions
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.7 , NEMO-consortium (2014)
43   !! $Id: sbcabl.F90 6416 2016-04-01 12:22:17Z clem $
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE sbc_abl_init
49      !!---------------------------------------------------------------------
50      !!                    ***  ROUTINE sbc_abl_init  ***
51      !!
52      !! ** Purposes :   - read namelist section namsbc_abl
53      !!                 - initialize and check parameter values
54      !!                 - initialize variables of ABL model
55      !!
56      !!----------------------------------------------------------------------
57      INTEGER            ::   ji, jj, jk, jbak, jbak_dta       ! dummy loop indices
58      INTEGER            ::   ios, ierror, ioptio              ! Local integer
59      INTEGER            ::   inum, indims, idimsz(4), id
60      CHARACTER(len=100) ::   cn_dir, cn_dom                   ! Atmospheric grid directory
61      REAL(wp)           ::   zcff,zcff1
62      LOGICAL            ::   lluldl
63      NAMELIST/namsbc_abl/ cn_dir , cn_dom, ln_hpgls_frc, ln_geos_winds,          &
64         &                 nn_dyn_restore,                                        &
65         &                 rn_ldyn_min , rn_ldyn_max, rn_ltra_min, rn_ltra_max,   &
66         &                 nn_amxl, rn_cm, rn_ct, rn_ce, rn_ceps, rn_Rod, rn_Ric, &
67         &                 ln_smth_pblh       
68      !!---------------------------------------------------------------------
69
70      REWIND( numnam_ref )              ! Namelist namsbc_abl in reference namelist : ABL parameters
71      READ  ( numnam_ref, namsbc_abl, IOSTAT = ios, ERR = 901 )
72901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in reference namelist', lwp )
73      !
74      REWIND( numnam_cfg )              ! Namelist namsbc_abl in configuration namelist : ABL parameters
75      READ  ( numnam_cfg, namsbc_abl, IOSTAT = ios, ERR = 902 )
76902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in configuration namelist', lwp )
77      !
78      IF(lwm) WRITE( numond, namsbc_abl )
79      !
80      ! Check ABL mixing length option     
81      IF( nn_amxl  < 0   .OR.  nn_amxl  > 2 )   &
82         &                 CALL ctl_stop( 'abl_init : bad flag, nn_amxl must be  0, 1 or 2 ' )         
83      !
84      ! Check ABL dyn restore option         
85      IF( nn_dyn_restore  < 0   .OR.  nn_dyn_restore  > 2 )   &
86         &                 CALL ctl_stop( 'abl_init : bad flag, nn_dyn_restore must be  0, 1 or 2 ' )           
87      !
88      !!---------------------------------------------------------------------
89      !! Control prints
90      !!---------------------------------------------------------------------
91      IF(lwp) THEN                              ! Control print (other namelist variable)
92         WRITE(numout,*)
93         WRITE(numout,*) '      ABL -- cn_dir         = ', cn_dir
94         WRITE(numout,*) '      ABL -- cn_dom         = ', cn_dom
95         IF( ln_hpgls_frc ) THEN
96            WRITE(numout,*) '      ABL -- winds forced by large-scale pressure gradient'
97            IF(ln_geos_winds) THEN
98               ln_geos_winds = .FALSE.
99               WRITE(numout,*) '      ABL -- geostrophic guide disabled (not compatible with ln_hpgls_frc = .T.)'                 
100            END IF
101         ELSE IF( ln_geos_winds ) THEN
102            WRITE(numout,*) '      ABL -- winds forced by geostrophic winds'               
103         ELSE
104            WRITE(numout,*) '      ABL -- Geostrophic winds and large-scale pressure gradient are ignored'               
105         END IF
106         !
107         SELECT CASE ( nn_dyn_restore )
108         CASE ( 0 )   
109            WRITE(numout,*) '      ABL -- No restoring for ABL winds'   
110         CASE ( 1 )
111            WRITE(numout,*) '      ABL -- Restoring of ABL winds only in the equatorial region '             
112         CASE ( 2 )
113            WRITE(numout,*) '      ABL -- Restoring of ABL winds activated everywhere '             
114         END SELECT
115         !
116       IF( ln_smth_pblh )  WRITE(numout,*) '      ABL -- Smoothing of PBL height is activated'   
117       !
118      ENDIF
119
120      !!---------------------------------------------------------------------
121      !! Convert nudging coefficient from hours to 1/sec
122      !!---------------------------------------------------------------------         
123      zcff        = 1._wp / 3600._wp
124      rn_ldyn_min = zcff / rn_ldyn_min
125      rn_ldyn_max = zcff / rn_ldyn_max     
126      rn_ltra_min = zcff / rn_ltra_min
127      rn_ltra_max = zcff / rn_ltra_max   
128
129      !!---------------------------------------------------------------------
130      !! ABL grid initialization
131      !!---------------------------------------------------------------------     
132      CALL iom_open( TRIM(cn_dir)//TRIM(cn_dom), inum ) 
133      id     = iom_varid( inum, 'e3t_abl', kdimsz=idimsz, kndims=indims, lduld=lluldl )
134      jpka   = idimsz(indims - COUNT( (/lluldl/) ) ) 
135      jpkam1 = jpka - 1
136
137      IF( abl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'abl_init : unable to allocate arrays' )
138      CALL iom_get( inum, jpdom_unknown, 'e3t_abl', e3t_abl(:) )
139      CALL iom_get( inum, jpdom_unknown, 'e3w_abl', e3w_abl(:) )
140      CALL iom_get( inum, jpdom_unknown, 'ght_abl', ght_abl(:) )
141      CALL iom_get( inum, jpdom_unknown, 'ghw_abl', ghw_abl(:) )
142      CALL iom_close( inum )
143
144      IF(lwp) THEN
145         WRITE(numout,*)
146         WRITE(numout,*) '    sbc_abl_init   : ABL Reference vertical grid'
147         WRITE(numout,*) '    ~~~~~~~'     
148         WRITE(numout, "(9x,'  level ght_abl   ghw_abl   e3t_abl   e3w_abl  ')" )
149         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, ght_abl(jk), ghw_abl(jk), e3t_abl(jk), e3w_abl(jk), jk = 1, jpka )
150      END IF
151
152      !!---------------------------------------------------------------------
153      !! Check TKE closure parameters
154      !!---------------------------------------------------------------------
155      rn_Sch  = rn_ce / rn_cm
156      mxl_min = (avm_bak / rn_cm) / sqrt( tke_min )
157
158      IF(lwp) THEN
159         WRITE(numout,*)
160         WRITE(numout,*) '    abl_zdf_tke   : ABL TKE turbulent closure'
161         WRITE(numout,*) '    ~~~~~~~~~~~'     
162         IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale '
163         IF(nn_amxl==1) WRITE(numout,*) 'length-scale based on the distance to the PBL height '
164         WRITE(numout,*) ' Minimum value of atmospheric TKE           = ',tke_min,' m^2 s^-2'
165         WRITE(numout,*) ' Minimum value of atmospheric mixing length = ',mxl_min,' m'
166         WRITE(numout,*) ' Constant for turbulent viscosity           = ',rn_Cm         
167         WRITE(numout,*) ' Constant for turbulent diffusivity         = ',rn_Ct         
168         WRITE(numout,*) ' Constant for Schmidt number                = ',rn_Sch           
169         WRITE(numout,*) ' Constant for TKE dissipation               = ',rn_Ceps                       
170      END IF
171
172      !!-------------------------------------------------------------------------------------------
173      !! Compute parameters to build the vertical profile for the nudging term (used in abl_stp())
174      !!-------------------------------------------------------------------------------------------       
175      zcff1       = 1._wp / ( jp_bmax - jp_bmin )**3
176      ! for active tracers
177      jp_alp3_tra = -2._wp * zcff1                       * ( rn_ltra_max - rn_ltra_min )
178      jp_alp2_tra =  3._wp * zcff1 * (jp_bmax + jp_bmin) * ( rn_ltra_max - rn_ltra_min )
179      jp_alp1_tra = -6._wp * zcff1 *  jp_bmax * jp_bmin  * ( rn_ltra_max - rn_ltra_min )
180      jp_alp0_tra = zcff1 * (  rn_ltra_max * jp_bmin*jp_bmin * (3._wp*jp_bmax - jp_bmin)      &
181         &                   - rn_ltra_min * jp_bmax*jp_bmax * (3._wp*jp_bmin - jp_bmax) ) 
182      ! for dynamics
183      jp_alp3_dyn = -2._wp * zcff1                       * ( rn_ldyn_max - rn_ldyn_min )
184      jp_alp2_dyn =  3._wp * zcff1 * (jp_bmax + jp_bmin) * ( rn_ldyn_max - rn_ldyn_min )
185      jp_alp1_dyn = -6._wp * zcff1 *  jp_bmax * jp_bmin  * ( rn_ldyn_max - rn_ldyn_min )
186      jp_alp0_dyn = zcff1 * (  rn_ldyn_max * jp_bmin*jp_bmin * (3._wp*jp_bmax - jp_bmin)      &
187         &                   - rn_ldyn_min * jp_bmax*jp_bmax * (3._wp*jp_bmin - jp_bmax) )       
188
189      jp_pblh_min = ghw_abl(     4) / jp_bmin  !<-- at least 3 grid points at the bottom have value rn_ltra_min
190      jp_pblh_max = ghw_abl(jpka-3) / jp_bmax  !<-- at least 3 grid points at the top    have value rn_ltra_max
191
192      ! Check parameters for dynamics
193      zcff  = jp_alp3_dyn * jp_bmin**3 + jp_alp2_dyn * jp_bmin**2   &
194         &  + jp_alp1_dyn * jp_bmin    + jp_alp0_dyn 
195      zcff1 = jp_alp3_dyn * jp_bmax**3 + jp_alp2_dyn * jp_bmax**2   &
196         &  + jp_alp1_dyn * jp_bmax    + jp_alp0_dyn     
197      IF(lwp) THEN
198         IF(nn_dyn_restore > 0) THEN
199            WRITE(numout,*) ' ABL Minimum value for dynamics restoring = ',zcff  * rdt
200            WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1 * rdt
201            ! Check that restoring coefficients are between 0 and 1
202            IF( zcff1 * rdt > 1._wp .OR. zcff1 * rdt < 0._wp )   &
203               &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_max' )
204            IF( zcff  * rdt > 1._wp .OR. zcff  * rdt < 0._wp )   &
205               &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_min' )
206            IF( zcff  * rdt > zcff1 * rdt ) &
207               &                   CALL ctl_stop( 'abl_init : rn_ldyn_max should be smaller than rn_ldyn_min' )
208         END IF
209      END IF
210
211      ! Check parameters for active tracers
212      zcff  = jp_alp3_tra * jp_bmin**3 + jp_alp2_tra * jp_bmin**2   &
213         &  + jp_alp1_tra * jp_bmin    + jp_alp0_tra 
214      zcff1 = jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2   &
215         &  + jp_alp1_tra * jp_bmax    + jp_alp0_tra 
216      IF(lwp) THEN
217         WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff  * rdt
218         WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1 * rdt
219         ! Check that restoring coefficients are between 0 and 1
220         IF( zcff1 * rdt > 1._wp .OR. zcff1 * rdt < 0._wp )   &
221            &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_max' )
222         IF( zcff  * rdt > 1._wp .OR. zcff  * rdt < 0._wp )   &
223            &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_min' )
224         IF( zcff  * rdt > zcff1  * rdt ) &
225            &                   CALL ctl_stop( 'abl_init : rn_ltra_max should be smaller than rn_ltra_min' )
226      END IF
227
228      !!-------------------------------------------------------------------------------------------
229      !! Initialize Coriolis frequency, equatorial restoring and land/sea mask
230      !!------------------------------------------------------------------------------------------- 
231      fft_abl(:,:) = 2._wp * omega * SIN( rad * gphit(:,:) )
232
233      ! Equatorial restoring
234      IF( nn_dyn_restore == 1 ) THEN
235         zcff         = 2._wp * omega * SIN( rad * 90._wp )   !++ fmax
236         rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff)/zcff ) )**8
237      ELSE
238         rest_eq(:,:) = 1._wp
239      END IF
240      ! T-mask
241      msk_abl(:,:) = tmask(:,:,1)
242
243      !!-------------------------------------------------------------------------------------------
244
245      ! initialize 2D bulk fields AND 3D abl data
246      CALL sbc_blk_init
247
248      ! initialize ABL from data or restart
249      !!GS  disabled for now
250      !IF( ln_rstart ) THEN
251      !  CALL ctl_stop( 'STOP', 'sbc_abl_init: restart mode not supported yet' )
252      !ELSE
253
254      CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step
255
256      ! Initialize the time index for now time (nt_n) and after time (nt_a)
257      nt_n = 1 + MOD( nit000  , 2)
258      nt_a = 1 + MOD( nit000+1, 2)
259
260       u_abl(:,:,:,nt_n      ) = sf(jp_wndi)%fnow(:,:,:)
261       v_abl(:,:,:,nt_n      ) = sf(jp_wndj)%fnow(:,:,:)
262      tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:)
263      tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:)
264
265      tke_abl(:,:,:,nt_n     ) = tke_min
266      avm_abl(:,:,:          ) = avm_bak
267      avt_abl(:,:,:          ) = avt_bak
268      mxl_abl(:,:,:          ) = mxl_min 
269      pblh   (:,:            ) = ghw_abl( 3 )  !<-- assume that the pbl contains 3 grid points
270      u_abl  (:,:,:,nt_a     ) = 0._wp
271      v_abl  (:,:,:,nt_a     ) = 0._wp
272      tq_abl (:,:,:,nt_a,:   ) = 0._wp
273      tke_abl(:,:,:,nt_a     ) = 0._wp
274      !ENDIF
275      !!GS restart case not supported
276     
277   END SUBROUTINE sbc_abl_init
278
279 
280
281   SUBROUTINE sbc_abl( kt )
282      !!---------------------------------------------------------------------
283      !!                     ***  ROUTINE sbc_abl  ***
284      !!
285      !! ** Purpose :   provide the momentum, heat and freshwater fluxes at
286      !!      the ocean surface from an ABL calculation at each oceanic time step 
287      !!
288      !! ** Method  :
289      !!              - Pre-compute part of turbulent fluxes in blk_oce_1
290      !!              - Perform 1 time-step of the ABL model   
291      !!              - Finalize flux computation in blk_oce_2
292      !!
293      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
294      !!              - vtau    : j-component of the stress at V-point  (N/m2)
295      !!              - taum    : Wind stress module at T-point         (N/m2)
296      !!              - wndm    : Wind speed module at T-point          (m/s)
297      !!              - qsr     : Solar heat flux over the ocean        (W/m2)
298      !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
299      !!              - emp     : evaporation minus precipitation       (kg/m2/s)
300      !!
301      !!---------------------------------------------------------------------
302      INTEGER ,         INTENT(in) ::   kt   ! ocean time step
303      !!
304      REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zevp
305      INTEGER                      ::   jbak, jbak_dta, ji, jj
306      !!---------------------------------------------------------------------
307      !
308      !!-------------------------------------------------------------------------------------------
309      !! 1 - Read Atmospheric 3D data for large-scale forcing
310      !!-------------------------------------------------------------------------------------------
311     
312      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
313     
314      !!-------------------------------------------------------------------------------------------
315      !! 2 - Compute Cd x ||U||, Ch x ||U||, Ce x ||U||, and SSQ using now fields
316      !!------------------------------------------------------------------------------------------- 
317         
318      CALL blk_oce_1( kt, u_abl(:,:,2,nt_n), v_abl(:,:,2,nt_n), tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),             &
319              &           sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, &
320              &              zssq, zcd_du, zsen, zevp )
321 
322      !!-------------------------------------------------------------------------------------------
323      !! 3 - Advance ABL variables from now (n) to after (n+1)
324      !!-------------------------------------------------------------------------------------------   
325 
326      CALL abl_stp( kt, sst_m, ssu_m, ssv_m, zssq, &                            ! in
327              &         sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:),   &
328              &         sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:), sf(jp_slp)%fnow(:,:,1),  &
329              &         sf(jp_hpgi)%fnow(:,:,:), sf(jp_hpgj)%fnow(:,:,:),                          &
330              &         zcd_du, zsen, zevp,    &                                ! in/out
331              &         wndm, utau, vtau, taum )                                ! out
332
333      !!-------------------------------------------------------------------------------------------
334      !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since
335      !!                                                                time swap is done in abl_stp
336      !!-------------------------------------------------------------------------------------------   
337
338      CALL blk_oce_2( kt, tq_abl(:,:,2,nt_n,jp_ta), sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1),     &
339              &           sf(jp_prec)%fnow(:,:,1),  sf(jp_snow)%fnow(:,:,1),   &
340              &           sst_m , zsen, zevp )
341
342   END SUBROUTINE sbc_abl
343
344   !!======================================================================
345END MODULE sbcabl
Note: See TracBrowser for help on using the repository browser.