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 NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/SBC – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/SBC/sbcflx.F90 @ 15675

Last change on this file since 15675 was 15675, checked in by hadjt, 2 years ago

Bugfix SBC/sbcflx.F90:

Enda emailed (27 January 2022 09:28) to say a bug had been spotted in SBC/sbcflx.F90, where the wind stress was halved.

I have removed this 0.5 factor in (my) line 250.
File size: 14.0 KB
Line 
1MODULE sbcflx
2   !!======================================================================
3   !!                       ***  MODULE  sbcflx  ***
4   !! Ocean forcing:  momentum, heat and freshwater flux formulation
5   !!=====================================================================
6   !! History :  1.0  !  2006-06  (G. Madec)  Original code
7   !!            3.3  !  2010-10  (S. Masson)  add diurnal cycle
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   namflx   : flux formulation namlist
12   !!   sbc_flx  : flux formulation as ocean surface boundary condition (forced mode, fluxes read in NetCDF files)
13   !!----------------------------------------------------------------------
14   USE oce             ! ocean dynamics and tracers
15   USE dom_oce         ! ocean space and time domain
16   USE sbc_oce         ! surface boundary condition: ocean fields
17   USE sbcdcy          ! surface boundary condition: diurnal cycle on qsr
18   USE phycst          ! physical constants
19   !
20   USE fldread         ! read input fields
21   USE iom             ! IOM library
22   USE in_out_manager  ! I/O manager
23   USE lib_mpp         ! distribued memory computing library
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC sbc_flx       ! routine called by step.F90
30
31   INTEGER , PARAMETER ::   jp_utau = 1   ! index of wind stress (i-component) file
32   INTEGER , PARAMETER ::   jp_vtau = 2   ! index of wind stress (j-component) file
33   INTEGER , PARAMETER ::   jp_qtot = 3   ! index of total (non solar+solar) heat file
34   INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file
35   INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file
36 !!INTEGER , PARAMETER ::   jp_sfx  = 6   ! index of salt flux flux
37   INTEGER , PARAMETER ::   jp_press = 6  ! index of pressure for UKMO shelf fluxes
38   INTEGER , PARAMETER ::   jpfld   = 6 !! 6 ! maximum number of files to read
39
40   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
41   ! JT
42   LOGICAL , PUBLIC    ::   ln_shelf_flx = .TRUE. ! UKMO SHELF specific flux flag
43   INTEGER             ::   jpfld_local   ! maximum number of files to read (locally modified depending on ln_shelf_flx)
44   ! JT
45
46
47   !! * Substitutions
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE sbc_flx( kt )
57      !!---------------------------------------------------------------------
58      !!                    ***  ROUTINE sbc_flx  ***
59      !!                   
60      !! ** Purpose :   provide at each time step the surface ocean fluxes
61      !!                (momentum, heat, freshwater and runoff)
62      !!
63      !! ** Method  : - READ each fluxes in NetCDF files:
64      !!                   i-component of the stress              utau  (N/m2)
65      !!                   j-component of the stress              vtau  (N/m2)
66      !!                   net downward heat flux                 qtot  (watt/m2)
67      !!                   net downward radiative flux            qsr   (watt/m2)
68      !!                   net upward freshwater (evapo - precip) emp   (kg/m2/s)
69      !!                   salt flux                              sfx   (pss*dh*rho/dt => g/m2/s)
70      !!
71      !!      CAUTION :  - never mask the surface stress fields
72      !!                 - the stress is assumed to be in the (i,j) mesh referential
73      !!
74      !! ** Action  :   update at each time-step
75      !!              - utau, vtau  i- and j-component of the wind stress
76      !!              - taum        wind stress module at T-point
77      !!              - wndm        10m wind module at T-point
78      !!              - qns         non solar heat flux including heat flux due to emp
79      !!              - qsr         solar heat flux
80      !!              - emp         upward mass flux (evap. - precip.)
81      !!              - sfx         salt flux; set to zero at nit000 but possibly non-zero
82      !!                            if ice
83      !!----------------------------------------------------------------------
84      INTEGER, INTENT(in) ::   kt   ! ocean time step
85      !!
86      INTEGER  ::   ji, jj, jf            ! dummy indices
87      INTEGER  ::   ierror                ! return error code
88      INTEGER  ::   ios                   ! Local integer output status for namelist read
89      REAL(wp) ::   zfact                 ! temporary scalar
90      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
91      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
92      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
93      ! JT
94      REAL     ::   cs                    ! UKMO SHELF: Friction co-efficient at surface
95      REAL     ::   totwindspd            ! UKMO SHELF: Magnitude of wind speed vector
96
97      REAL(wp) ::   rhoa  = 1.22         ! Air density kg/m3
98      REAL(wp) ::   cdrag = 1.5e-3       ! drag coefficient
99      ! JT
100      !!
101      CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files
102      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                    ! array of namelist information structures
103      ! JT
104      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp, sn_press !!, sn_sfx ! informations about the fields to be read
105      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp,   &
106      &                    sn_press, ln_shelf_flx
107      ! JT
108      !JT remove nul from output.namelist.dyn
109      sn_utau%lname = ''
110      sn_vtau%lname = ''
111      sn_qtot%lname = ''
112      sn_qsr%lname = ''
113      sn_emp%lname = ''
114      sn_press%lname = ''
115      !JT remove nul from output.namelist.dyn
116      !!---------------------------------------------------------------------
117      !
118      IF( kt == nit000 ) THEN                ! First call kt=nit000 
119         ! set file information
120         REWIND( numnam_ref )              ! Namelist namsbc_flx in reference namelist : Files for fluxes
121         READ  ( numnam_ref, namsbc_flx, IOSTAT = ios, ERR = 901)
122901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_flx in reference namelist' )
123
124         REWIND( numnam_cfg )              ! Namelist namsbc_flx in configuration namelist : Files for fluxes
125         READ  ( numnam_cfg, namsbc_flx, IOSTAT = ios, ERR = 902 )
126902      IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_flx in configuration namelist' )
127         IF(lwm) WRITE ( numond, namsbc_flx ) 
128         
129         
130         !
131         !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing?
132         IF( ln_dm2dc .AND. sn_qsr%freqh /= 24. )   &
133            &   CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' ) 
134         !
135         !                                         ! store namelist information in an array
136         slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau
137         slf_i(jp_qtot) = sn_qtot   ;   slf_i(jp_qsr ) = sn_qsr 
138         slf_i(jp_emp ) = sn_emp !! ;   slf_i(jp_sfx ) = sn_sfx
139         !
140            ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
141            IF( ln_shelf_flx ) slf_i(jp_press) = sn_press
142   
143            ! define local jpfld depending on shelf_flx logical
144            IF( ln_shelf_flx ) THEN
145               jpfld_local = jpfld
146            ELSE
147               jpfld_local = jpfld-1
148            ENDIF
149            !
150         IF( ierror > 0 ) THEN   
151            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN 
152         ENDIF
153         DO ji= 1, jpfld_local
154            ALLOCATE( sf(ji)%fnow(jpi,jpj,1) )
155            IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) )
156         END DO
157         !                                         ! fill sf with slf_i and control print
158         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
159         !
160         sfx(:,:) = 0.0_wp                         ! salt flux due to freezing/melting (non-zero only if ice is present; set in limsbc(_2).F90)
161         !
162      ENDIF
163
164      CALL fld_read( kt, nn_fsbc, sf )                            ! input fields provided at the current time-step
165     
166      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency
167
168         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)  ! modify now Qsr to include the diurnal cycle
169         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
170         ENDIF
171
172            !!UKMO SHELF effect of atmospheric pressure on SSH
173            ! If using ln_apr_dyn, this is done there so don't repeat here.
174            !IF( ln_shelf_flx .AND. .NOT. ln_apr_dyn) THEN
175            !   DO jj = 1, jpjm1
176            !      DO ji = 1, jpim1
177            !         apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj)
178            !         apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj)
179            !      END DO
180            !   END DO
181            !ENDIF
182
183         DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
184            DO ji = 1, jpi
185               !JT
186               IF( ln_shelf_flx ) THEN
187                  !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing
188                  pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1)
189                  !! UKMO SHELF flux files contain wind speed not wind stress
190                  totwindspd = sqrt((sf(jp_utau)%fnow(ji,jj,1))**2.0 + (sf(jp_vtau)%fnow(ji,jj,1))**2.0)
191                  cs = 0.63 + (0.066 * totwindspd)
192                  utau(ji,jj) = (cs * (rhoa/rau0) * sf(jp_utau)%fnow(ji,jj,1) * totwindspd)* umask(ji,jj,1)
193                  vtau(ji,jj) = (cs * (rhoa/rau0) * sf(jp_vtau)%fnow(ji,jj,1) * totwindspd)* vmask(ji,jj,1)
194               ELSE
195                  utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)* umask(ji,jj,1)
196                  vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)* vmask(ji,jj,1)
197               ENDIF
198               !JT                   
199               qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1)                                * tmask(ji,jj,1)
200               IF( ln_shelf_flx ) THEN
201                  !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot)
202                  qns (ji,jj) = (sf(jp_qtot)%fnow(ji,jj,1)) * tmask(ji,jj,1)
203                  !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P
204                  emp (ji,jj) = (-1. * sf(jp_emp )%fnow(ji,jj,1)) * tmask(ji,jj,1)
205               ELSE
206                  qns (ji,jj) = (sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)) * tmask(ji,jj,1)
207                  emp (ji,jj) = (sf(jp_emp )%fnow(ji,jj,1)) * tmask(ji,jj,1)
208               ENDIF
209!               !utau(ji,jj) =   sf(jp_utau)%fnow(ji,jj,1)                              * umask(ji,jj,1)
210!               !vtau(ji,jj) =   sf(jp_vtau)%fnow(ji,jj,1)                              * vmask(ji,jj,1)
211!               !JT qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
212!               qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - qsr (ji,jj)              ) * tmask(ji,jj,1)
213!               emp (ji,jj) =   sf(jp_emp )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
214!               !!sfx (ji,jj) = sf(jp_sfx )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
215            END DO
216         END DO
217         !                                                        ! add to qns the heat due to e-p
218         !clem: I do not think it is needed
219         !!qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
220         !
221         ! clem: without these lbc calls, it seems that the northfold is not ok (true in 3.6, not sure in 4.x)
222
223         !JT Enabling as used in previous version using these climate forcings
224         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
225         !JT
226
227         CALL lbc_lnk_multi( 'sbcflx', utau, 'U', -1._wp, vtau, 'V', -1._wp, &
228            &                           qns, 'T',  1._wp, emp , 'T',  1._wp, qsr, 'T', 1._wp ) !! sfx, 'T', 1._wp  )
229         !
230         IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
231            WRITE(numout,*) 
232            WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK'
233            DO jf = 1, jpfld_local
234               IF( jf == jp_utau .OR. jf == jp_vtau )   zfact =     1.
235               IF( jf == jp_qtot .OR. jf == jp_qsr  )   zfact =     0.1
236               IF( jf == jp_emp                     )   zfact = 86400.
237               WRITE(numout,*) 
238               WRITE(numout,*) ' day: ', ndastp , TRIM(sf(jf)%clvar), ' * ', zfact
239            END DO
240         ENDIF
241         !
242      ENDIF
243      !                                                           ! module of wind stress and wind speed at T-point
244      ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
245      zcoef = 1. / ( zrhoa * zcdrag )
246      DO jj = 2, jpjm1
247         DO ji = fs_2, fs_jpim1   ! vect. opt.
248            ztx = ( utau(ji-1,jj  ) + utau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( umask(ji-1,jj  ,1), umask(ji,jj,1) ) )
249            zty = ( vtau(ji  ,jj-1) + vtau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( vmask(ji  ,jj-1,1), vmask(ji,jj,1) ) ) 
250            zmod = SQRT( ztx * ztx + zty * zty ) * tmask(ji,jj,1)
251            taum(ji,jj) = zmod
252            wndm(ji,jj) = SQRT( zmod * zcoef )  !!clem: not used?
253         END DO
254      END DO
255      !
256      CALL lbc_lnk_multi( 'sbcflx', taum, 'T', 1._wp, wndm, 'T', 1._wp )
257      !
258      !
259   END SUBROUTINE sbc_flx
260
261   !!======================================================================
262END MODULE sbcflx
Note: See TracBrowser for help on using the repository browser.