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.
sbcflx.F90 in branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90

Last change on this file was 12654, checked in by jcastill, 4 years ago

Further change for Met Office utils ticket 334: assume that the wind forcing components are in the U,V grid instead of the T grid

File size: 16.3 KB
RevLine 
[888]1MODULE sbcflx
2   !!======================================================================
3   !!                       ***  MODULE  sbcflx  ***
4   !! Ocean forcing:  momentum, heat and freshwater flux formulation
5   !!=====================================================================
[2528]6   !! History :  1.0  !  2006-06  (G. Madec)  Original code
7   !!            3.3  !  2010-10  (S. Masson)  add diurnal cycle
[888]8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   namflx   : flux formulation namlist
[2528]12   !!   sbc_flx  : flux formulation as ocean surface boundary condition (forced mode, fluxes read in NetCDF files)
[888]13   !!----------------------------------------------------------------------
14   USE oce             ! ocean dynamics and tracers
15   USE dom_oce         ! ocean space and time domain
[2528]16   USE sbc_oce         ! surface boundary condition: ocean fields
17   USE sbcdcy          ! surface boundary condition: diurnal cycle on qsr
[888]18   USE phycst          ! physical constants
19   USE fldread         ! read input fields
20   USE iom             ! IOM library
21   USE in_out_manager  ! I/O manager
[11277]22   USE sbcwave         ! wave physics
[888]23   USE lib_mpp         ! distribued memory computing library
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[8416]25   USE wrk_nemo        ! work arrays
[888]26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC sbc_flx       ! routine called by step.F90
31
[8059]32   INTEGER , PARAMETER ::   jpfld   = 6   ! maximum number of files to read
[888]33   INTEGER , PARAMETER ::   jp_utau = 1   ! index of wind stress (i-component) file
34   INTEGER , PARAMETER ::   jp_vtau = 2   ! index of wind stress (j-component) file
35   INTEGER , PARAMETER ::   jp_qtot = 3   ! index of total (non solar+solar) heat file
36   INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file
37   INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file
[8059]38   INTEGER , PARAMETER ::   jp_press = 6  ! index of pressure for UKMO shelf fluxes
[888]39   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
[8059]40   LOGICAL , PUBLIC    ::   ln_shelf_flx = .FALSE. ! UKMO SHELF specific flux flag
[8416]41   LOGICAL , PUBLIC    ::   ln_rel_wind = .FALSE.  ! UKMO SHELF specific flux flag - relative winds
42   REAL(wp)            ::   rn_wfac                ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
[8059]43   INTEGER             ::   jpfld_local   ! maximum number of files to read (locally modified depending on ln_shelf_flx)
[888]44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
[1029]47#  include "vectopt_loop_substitute.h90"
[888]48   !!----------------------------------------------------------------------
[2528]49   !! NEMO/OPA 3.3 , NEMO-consortium (2010)
[1156]50   !! $Id$
[2715]51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[888]52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE sbc_flx( kt )
56      !!---------------------------------------------------------------------
57      !!                    ***  ROUTINE sbc_flx  ***
58      !!                   
59      !! ** Purpose :   provide at each time step the surface ocean fluxes
60      !!                (momentum, heat, freshwater and runoff)
61      !!
62      !! ** Method  : - READ each fluxes in NetCDF files:
63      !!                   i-component of the stress              utau  (N/m2)
64      !!                   j-component of the stress              vtau  (N/m2)
65      !!                   net downward heat flux                 qtot  (watt/m2)
66      !!                   net downward radiative flux            qsr   (watt/m2)
67      !!                   net upward freshwater (evapo - precip) emp   (kg/m2/s)
68      !!
69      !!      CAUTION :  - never mask the surface stress fields
[3625]70      !!                 - the stress is assumed to be in the (i,j) mesh referential
[888]71      !!
72      !! ** Action  :   update at each time-step
[1695]73      !!              - utau, vtau  i- and j-component of the wind stress
74      !!              - taum        wind stress module at T-point
75      !!              - wndm        10m wind module at T-point
[3625]76      !!              - qns         non solar heat flux including heat flux due to emp
77      !!              - qsr         solar heat flux
78      !!              - emp         upward mass flux (evap. - precip.)
79      !!              - sfx         salt flux; set to zero at nit000 but possibly non-zero
80      !!                            if ice is present (computed in limsbc(_2).F90)
[888]81      !!----------------------------------------------------------------------
82      INTEGER, INTENT(in) ::   kt   ! ocean time step
83      !!
[1695]84      INTEGER  ::   ji, jj, jf            ! dummy indices
85      INTEGER  ::   ierror                ! return error code
[4147]86      INTEGER  ::   ios                   ! Local integer output status for namelist read
[1695]87      REAL(wp) ::   zfact                 ! temporary scalar
88      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
89      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
[11277]90      REAL(wp) ::   totwind               ! UKMO SHELF: Module of wind speed
[1695]91      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
[8059]92      REAL     ::   cs                    ! UKMO SHELF: Friction co-efficient at surface
93      REAL     ::   totwindspd            ! UKMO SHELF: Magnitude of wind speed vector
[12654]94      REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at U,V-points
[8059]95
96      REAL(wp) ::   rhoa  = 1.22         ! Air density kg/m3
97      REAL(wp) ::   cdrag = 1.5e-3       ! drag coefficient
[888]98      !!
99      CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files
100      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                    ! array of namelist information structures
[8416]101      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp, sn_press  !  informations about the fields to be read
102      LOGICAL     ::   ln_foam_flx  = .FALSE.                     ! UKMO FOAM specific flux flag
103      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp,   & 
104      &                    ln_foam_flx, sn_press, ln_shelf_flx, ln_rel_wind,    &
105      &                    rn_wfac
[888]106      !!---------------------------------------------------------------------
[2528]107      !
108      IF( kt == nit000 ) THEN                ! First call kt=nit000 
[888]109         ! set file information
[4147]110         REWIND( numnam_ref )              ! Namelist namsbc_flx in reference namelist : Files for fluxes
111         READ  ( numnam_ref, namsbc_flx, IOSTAT = ios, ERR = 901)
112901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_flx in reference namelist', lwp )
113
114         REWIND( numnam_cfg )              ! Namelist namsbc_flx in configuration namelist : Files for fluxes
115         READ  ( numnam_cfg, namsbc_flx, IOSTAT = ios, ERR = 902 )
116902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_flx in configuration namelist', lwp )
[4624]117         IF(lwm) WRITE ( numond, namsbc_flx ) 
[2528]118         !
[11277]119         IF(lwp) THEN                        ! Namelist print
120            WRITE(numout,*) 
121            WRITE(numout,*) 'sbc_flx : Flux forcing' 
122            WRITE(numout,*) '~~~~~~~~~~~' 
123            WRITE(numout,*) '       Namelist namsbc_flx : shelf seas configuration (force with winds instead of momentum)' 
124            WRITE(numout,*) '          shelf seas configuration    ln_shelf_flx    = ', ln_shelf_flx 
125            WRITE(numout,*) '          relative wind speed         ln_rel_wind     = ', ln_rel_wind 
126            WRITE(numout,*) '          wind multiplication factor  rn_wfac         = ', rn_wfac 
127         ENDIF
[2528]128         !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing?
129         IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 )   &
130            &   CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' ) 
131         !
132         !                                         ! store namelist information in an array
[888]133         slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau
134         slf_i(jp_qtot) = sn_qtot   ;   slf_i(jp_qsr ) = sn_qsr 
135         slf_i(jp_emp ) = sn_emp
[2528]136         !
[12599]137         ! define local jpfld depending on shelf_flx logical
138         IF( ln_shelf_flx ) THEN
139            slf_i(jp_press) = sn_press
140            jpfld_local = jpfld
141         ELSE
142            jpfld_local = jpfld-1
143         ENDIF
144         !
145         ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
146         !
[1133]147         IF( ierror > 0 ) THEN   
148            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN 
[888]149         ENDIF
[8059]150         DO ji= 1, jpfld_local
[2528]151            ALLOCATE( sf(ji)%fnow(jpi,jpj,1) )
152            IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) )
[1200]153         END DO
[2528]154         !                                         ! fill sf with slf_i and control print
[1133]155         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
[888]156         !
[3625]157         sfx(:,:) = 0.0_wp                         ! salt flux due to freezing/melting (non-zero only if ice is present; set in limsbc(_2).F90)
158         !
[888]159      ENDIF
160
[2528]161      CALL fld_read( kt, nn_fsbc, sf )                            ! input fields provided at the current time-step
162     
163      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency
[888]164
[12599]165         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )   ! modify now Qsr to include the diurnal cycle
166         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)
167         ENDIF
168
[12654]169         !! UKMO SHELF flux files contain wind speed not wind stress
[8416]170         IF( ln_shelf_flx ) THEN
171            CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j )
172
[12599]173            ! set the ocean fluxes from read fields
174            DO jj = 1, jpj
175               DO ji = 1, jpi
176                  !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing
177                  pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1)
[12654]178
179                  qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1)
180                  !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot)
181                  qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1)
182                  !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P
183                  emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1)
184
185                  !!UKMO SHELF wind speed relative to surface currents - put here to allow merging with coupling branch
[12599]186                  IF( ln_rel_wind ) THEN
[8918]187                     zwnd_i(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) - rn_wfac * ssu_m(ji,jj)
188                     zwnd_j(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) - rn_wfac * ssv_m(ji,jj)
[12599]189                  ELSE
190                     zwnd_i(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
191                     zwnd_j(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
192                  ENDIF
[12654]193               END DO
194            END DO
[12599]195
[12654]196            ! Calculate wind speed from the wind components
197            DO jj = 2, jpjm1
198               DO ji = fs_2, fs_jpim1   ! vect. opt.
199                  ztx = zwnd_i(ji-1,jj  ) + zwnd_i(ji,jj) 
200                  zty = zwnd_j(ji  ,jj-1) + zwnd_j(ji,jj) 
201                  wndm(ji,jj) = 0.5 * SQRT( ztx * ztx + zty * zty )
202               END DO
203            END DO
204
205            ! add modification due to drag coefficient read from wave forcing
206            DO jj = 1, jpjm1
207               DO ji = 1, jpim1
[12599]208                  IF( ln_cdgw .AND. nn_drag == jp_std ) THEN
209                     IF( cpl_wdrag ) THEN 
210                        ! reset utau and vtau to the wind components: the momentum will
211                        ! be calculated from the coupled value of the drag coefficient
212                        utau(ji,jj) = zwnd_i(ji,jj) 
213                        vtau(ji,jj) = zwnd_j(ji,jj) 
214                     ELSE
[12654]215                        utau(ji,jj) = zrhoa * 0.25 * ( cdn_wave(ji,jj) + cdn_wave(ji+1,jj) ) * zwnd_i(ji,jj) * &
216                                                     ( wndm(ji,jj)     + wndm(ji+1,jj) )
217                        vtau(ji,jj) = zrhoa * 0.25 * ( cdn_wave(ji,jj) + cdn_wave(ji,jj+1) ) * zwnd_j(ji,jj) * &
218                                                     ( wndm(ji,jj)     + wndm(ji,jj+1) )
[12599]219                     ENDIF
220                  ELSE IF( nn_drag == jp_const ) THEN
[12654]221                     utau(ji,jj) = zrhoa * zcdrag * zwnd_i(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) )
222                     vtau(ji,jj) = zrhoa * zcdrag * zwnd_j(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) )
[12599]223                  ELSE IF( nn_drag == jp_ukmo ) THEN
[12654]224                     totwindspd = sqrt(zwnd_i(ji,jj)*zwnd_i(ji,jj) + zwnd_j(ji,jj)*zwnd_j(ji,jj))
225                     cs = 0.63 + (0.066 * totwindspd)
226                     utau(ji,jj) = cs * (rhoa/rau0) * zwnd_i(ji,jj) * totwindspd
227                     vtau(ji,jj) = cs * (rhoa/rau0) * zwnd_j(ji,jj) * totwindspd
[12599]228                  ENDIF 
229
230                  !!UKMO SHELF effect of atmospheric pressure on SSH
231                  ! If using ln_apr_dyn, this is done there so don't repeat here.
232                  IF( .NOT. ln_apr_dyn) THEN
[8059]233                     apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj)
234                     apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj)
[12599]235                  ENDIF
[8059]236               END DO
[1274]237            END DO
[12654]238
239            CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j )
[12599]240         ELSE
241            DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
242               DO ji = 1, jpi
243                  utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
244                  vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
245                  qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1)
246                  IF( ln_foam_flx ) THEN
247                     !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot)
248                     qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1)
249                     !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P
250                     emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1)
251                  ELSE
252                     qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)
253                     emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1)
254                  ENDIF
255               END DO
256            END DO
[12654]257         ENDIF
[12599]258            !                                                        ! module of wind stress and wind speed at T-point
[12654]259         zcoef = 1. / ( zrhoa * zcdrag )
[12599]260!CDIR NOVERRCHK
[12654]261         DO jj = 2, jpjm1
[12599]262!CDIR NOVERRCHK
[12654]263            DO ji = fs_2, fs_jpim1   ! vect. opt.
264               ztx = utau(ji-1,jj  ) + utau(ji,jj) 
265               zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
266               zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
267               taum(ji,jj) = zmod
268               IF( .NOT. ln_shelf_flx ) wndm(ji,jj) = SQRT( zmod * zcoef )
[12599]269            END DO
[12654]270         END DO
[3625]271         !                                                        ! add to qns the heat due to e-p
272         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
273         !
[12599]274         !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe
275         IF( ln_foam_flx .OR. ln_shelf_flx ) THEN
276            CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. ) 
277         ENDIF
[8059]278   
[4990]279         taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1)
[12599]280         CALL lbc_lnk_multi( taum(:,:), 'T', 1., wndm(:,:), 'T', 1. )
[1695]281
[2528]282         IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
[888]283            WRITE(numout,*) 
[1274]284            WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK'
[8059]285            DO jf = 1, jpfld_local
[1274]286               IF( jf == jp_utau .OR. jf == jp_vtau )   zfact =     1.
287               IF( jf == jp_qtot .OR. jf == jp_qsr  )   zfact =     0.1
288               IF( jf == jp_emp                     )   zfact = 86400.
289               WRITE(numout,*) 
290               WRITE(numout,*) ' day: ', ndastp , TRIM(sf(jf)%clvar), ' * ', zfact
291               CALL prihre( sf(jf)%fnow, jpi, jpj, 1, jpi, 20, 1, jpj, 10, zfact, numout )
292            END DO
293            CALL FLUSH(numout)
294         ENDIF
295         !
[888]296      ENDIF
297      !
298   END SUBROUTINE sbc_flx
299
300   !!======================================================================
301END MODULE sbcflx
Note: See TracBrowser for help on using the repository browser.