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.
sbcabl.F90 in NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL – NEMO

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

Last change on this file since 11858 was 11858, checked in by gsamson, 5 years ago

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