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 @ 15310

Last change on this file since 15310 was 15310, checked in by hadjt, 3 years ago

Wind velocities rather than wind stress

SBC/sbcflx.F90 edited to allow wind velocities rather than wind stresses.
namelist keyword added to namsbc_flx: sn_press, ln_shelf_flx

The version is committed in restrospect, so may not include all the changes to work.
File size: 13.1 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
109      !!---------------------------------------------------------------------
110      !
111      IF( kt == nit000 ) THEN                ! First call kt=nit000 
112         ! set file information
113         REWIND( numnam_ref )              ! Namelist namsbc_flx in reference namelist : Files for fluxes
114         READ  ( numnam_ref, namsbc_flx, IOSTAT = ios, ERR = 901)
115901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_flx in reference namelist' )
116
117         REWIND( numnam_cfg )              ! Namelist namsbc_flx in configuration namelist : Files for fluxes
118         READ  ( numnam_cfg, namsbc_flx, IOSTAT = ios, ERR = 902 )
119902      IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_flx in configuration namelist' )
120         IF(lwm) WRITE ( numond, namsbc_flx ) 
121         
122         
123         !
124         !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing?
125         IF( ln_dm2dc .AND. sn_qsr%freqh /= 24. )   &
126            &   CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' ) 
127         !
128         !                                         ! store namelist information in an array
129         slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau
130         slf_i(jp_qtot) = sn_qtot   ;   slf_i(jp_qsr ) = sn_qsr 
131         slf_i(jp_emp ) = sn_emp !! ;   slf_i(jp_sfx ) = sn_sfx
132         !
133            ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
134            IF( ln_shelf_flx ) slf_i(jp_press) = sn_press
135   
136            ! define local jpfld depending on shelf_flx logical
137            IF( ln_shelf_flx ) THEN
138               jpfld_local = jpfld
139            ELSE
140               jpfld_local = jpfld-1
141            ENDIF
142            !
143         IF( ierror > 0 ) THEN   
144            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN 
145         ENDIF
146         DO ji= 1, jpfld_local
147            ALLOCATE( sf(ji)%fnow(jpi,jpj,1) )
148            IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) )
149         END DO
150         !                                         ! fill sf with slf_i and control print
151         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
152         !
153         sfx(:,:) = 0.0_wp                         ! salt flux due to freezing/melting (non-zero only if ice is present; set in limsbc(_2).F90)
154         !
155      ENDIF
156
157      CALL fld_read( kt, nn_fsbc, sf )                            ! input fields provided at the current time-step
158     
159      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency
160
161         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)  ! modify now Qsr to include the diurnal cycle
162         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
163         ENDIF
164
165            !!UKMO SHELF effect of atmospheric pressure on SSH
166            ! If using ln_apr_dyn, this is done there so don't repeat here.
167            !IF( ln_shelf_flx .AND. .NOT. ln_apr_dyn) THEN
168            !   DO jj = 1, jpjm1
169            !      DO ji = 1, jpim1
170            !         apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj)
171            !         apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj)
172            !      END DO
173            !   END DO
174            !ENDIF
175
176         DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
177            DO ji = 1, jpi
178               !JT
179               IF( ln_shelf_flx ) THEN
180                  !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing
181                  pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1)
182                  !! UKMO SHELF flux files contain wind speed not wind stress
183                  totwindspd = sqrt((sf(jp_utau)%fnow(ji,jj,1))**2.0 + (sf(jp_vtau)%fnow(ji,jj,1))**2.0)
184                  cs = 0.63 + (0.066 * totwindspd)
185                  utau(ji,jj) = (cs * (rhoa/rau0) * sf(jp_utau)%fnow(ji,jj,1) * totwindspd)* umask(ji,jj,1)
186                  vtau(ji,jj) = (cs * (rhoa/rau0) * sf(jp_vtau)%fnow(ji,jj,1) * totwindspd)* vmask(ji,jj,1)
187               ELSE
188                  utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)* umask(ji,jj,1)
189                  vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)* vmask(ji,jj,1)
190               ENDIF
191               !JT                   
192               qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1)                                * tmask(ji,jj,1)
193               !utau(ji,jj) =   sf(jp_utau)%fnow(ji,jj,1)                              * umask(ji,jj,1)
194               !vtau(ji,jj) =   sf(jp_vtau)%fnow(ji,jj,1)                              * vmask(ji,jj,1)
195               !JT qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
196               qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - qsr (ji,jj)              ) * tmask(ji,jj,1)
197               emp (ji,jj) =   sf(jp_emp )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
198               !!sfx (ji,jj) = sf(jp_sfx )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
199            END DO
200         END DO
201         !                                                        ! add to qns the heat due to e-p
202         !clem: I do not think it is needed
203         !!qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
204         !
205         ! clem: without these lbc calls, it seems that the northfold is not ok (true in 3.6, not sure in 4.x)
206
207         !JT Enabling as used in previous version using these climate forcings
208         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
209         !JT
210
211         CALL lbc_lnk_multi( 'sbcflx', utau, 'U', -1._wp, vtau, 'V', -1._wp, &
212            &                           qns, 'T',  1._wp, emp , 'T',  1._wp, qsr, 'T', 1._wp ) !! sfx, 'T', 1._wp  )
213         !
214         IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
215            WRITE(numout,*) 
216            WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK'
217            DO jf = 1, jpfld_local
218               IF( jf == jp_utau .OR. jf == jp_vtau )   zfact =     1.
219               IF( jf == jp_qtot .OR. jf == jp_qsr  )   zfact =     0.1
220               IF( jf == jp_emp                     )   zfact = 86400.
221               WRITE(numout,*) 
222               WRITE(numout,*) ' day: ', ndastp , TRIM(sf(jf)%clvar), ' * ', zfact
223            END DO
224         ENDIF
225         !
226      ENDIF
227      !                                                           ! module of wind stress and wind speed at T-point
228      ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
229      zcoef = 1. / ( zrhoa * zcdrag )
230      DO jj = 2, jpjm1
231         DO ji = fs_2, fs_jpim1   ! vect. opt.
232            ztx = ( utau(ji-1,jj  ) + utau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( umask(ji-1,jj  ,1), umask(ji,jj,1) ) )
233            zty = ( vtau(ji  ,jj-1) + vtau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( vmask(ji  ,jj-1,1), vmask(ji,jj,1) ) ) 
234            zmod = 0.5_wp * SQRT( ztx * ztx + zty * zty ) * tmask(ji,jj,1)
235            taum(ji,jj) = zmod
236            wndm(ji,jj) = SQRT( zmod * zcoef )  !!clem: not used?
237         END DO
238      END DO
239      !
240      CALL lbc_lnk_multi( 'sbcflx', taum, 'T', 1._wp, wndm, 'T', 1._wp )
241      !
242      !
243   END SUBROUTINE sbc_flx
244
245   !!======================================================================
246END MODULE sbcflx
Note: See TracBrowser for help on using the repository browser.