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/r14075_ukmo_shelf/src/OCE/SBC – NEMO

source: NEMO/branches/UKMO/r14075_ukmo_shelf/src/OCE/SBC/sbcflx.F90 @ 15455

Last change on this file since 15455 was 15455, checked in by jcastill, 3 years ago

Code for uncoupled configurations, some changes for coupling may be needed yet - merged branch branches/UKMO/r14075_cpl-pressure@15423
File size: 14.5 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   ! maximum number of files to read
39
40   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
41   LOGICAL , PUBLIC    ::   ln_shelf_flx = .FALSE. ! UKMO SHELF specific flux flag
42   LOGICAL , PUBLIC    ::   ln_rel_wind  = .FALSE. ! UKMO SHELF specific flux flag - relative winds
43   REAL(wp)            ::   rn_wfac                ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
44   INTEGER             ::   jpfld_local   ! maximum number of files to read (locally modified depending on ln_shelf_flx)
45
46   !! * Substitutions
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
50   !! $Id$
51   !! Software governed by the CeCILL license (see ./LICENSE)
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      !!                   salt flux                              sfx   (pss*dh*rho/dt => g/m2/s)
69      !!
70      !!      CAUTION :  - never mask the surface stress fields
71      !!                 - the stress is assumed to be in the (i,j) mesh referential
72      !!
73      !! ** Action  :   update at each time-step
74      !!              - utau, vtau  i- and j-component of the wind stress
75      !!              - taum        wind stress module at T-point
76      !!              - wndm        10m wind module at T-point
77      !!              - qns         non solar heat flux including heat flux due to emp
78      !!              - qsr         solar heat flux
79      !!              - emp         upward mass flux (evap. - precip.)
80      !!              - sfx         salt flux; set to zero at nit000 but possibly non-zero
81      !!                            if ice
82      !!----------------------------------------------------------------------
83      INTEGER, INTENT(in) ::   kt   ! ocean time step
84      !!
85      INTEGER  ::   ji, jj, jf            ! dummy indices
86      INTEGER  ::   ierror                ! return error code
87      INTEGER  ::   ios                   ! Local integer output status for namelist read
88      REAL(wp) ::   zfact                 ! temporary scalar
89      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
90      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
91      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
92      REAL     ::   cs                    ! UKMO SHELF: Friction co-efficient at surface
93      REAL     ::   totwindspd            ! UKMO SHELF: Magnitude of wind speed vector
94      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
95   
96      REAL(wp) ::   rhoa  = 1.22         ! Air density kg/m3
97      REAL(wp) ::   cdrag = 1.5e-3       ! drag coefficient
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
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
106      !!---------------------------------------------------------------------
107      !
108      IF( kt == nit000 ) THEN                ! First call kt=nit000 
109         ! set file information
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' )
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' )
117         IF(lwm) WRITE ( numond, namsbc_flx ) 
118         !
119         !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing?
120         IF( ln_dm2dc .AND. sn_qsr%freqh /= 24. )   &
121            &   CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' ) 
122         !
123         !                                         ! store namelist information in an array
124         slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau
125         slf_i(jp_qtot) = sn_qtot   ;   slf_i(jp_qsr ) = sn_qsr 
126         slf_i(jp_emp ) = sn_emp !! ;   slf_i(jp_sfx ) = sn_sfx
127         !
128         IF( ln_shelf_flx ) slf_i(jp_press) = sn_press 
129 
130         ! define local jpfld depending on shelf_flx logical
131         IF( ln_shelf_flx ) THEN
132            jpfld_local = jpfld 
133         ELSE
134            jpfld_local = jpfld-1 
135         ENDIF 
136         !
137         ALLOCATE( sf(jpfld_local), STAT=ierror )        ! set sf structure
138         IF( ierror > 0 ) THEN   
139            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN 
140         ENDIF
141         DO ji= 1, jpfld_local
142            ALLOCATE( sf(ji)%fnow(jpi,jpj,1) )
143            IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) )
144         END DO
145         !                                         ! fill sf with slf_i and control print
146         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
147         !
148      ENDIF
149
150      CALL fld_read( kt, nn_fsbc, sf )                            ! input fields provided at the current time-step
151     
152      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency
153
154         !!UKMO SHELF wind speed relative to surface currents
155         IF( ln_shelf_flx ) THEN
156            ALLOCATE( zwnd_i(jpi,jpj), zwnd_j(jpi,jpj) )
157
158            IF( ln_rel_wind ) THEN
159               DO jj = 1, jpj
160                  DO ji = 1, jpi
161                     zwnd_i(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) - rn_wfac * ssu_m(ji,jj)
162                     zwnd_j(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) - rn_wfac * ssv_m(ji,jj)
163                  END DO
164               END DO
165            ELSE
166               zwnd_i(:,:) = sf(jp_utau)%fnow(:,:,1)
167               zwnd_j(:,:) = sf(jp_vtau)%fnow(:,:,1)
168            ENDIF
169         ENDIF
170
171         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)  ! modify now Qsr to include the diurnal cycle
172         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
173         ENDIF
174         !!UKMO SHELF effect of atmospheric pressure on SSH
175         ! If using ln_apr_dyn, this is done there so don't repeat here.
176         IF( ln_shelf_flx .AND. .NOT. ln_apr_dyn) THEN
177            DO jj = 1, jpjm1 
178               DO ji = 1, jpim1 
179                  apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj) 
180                  apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj) 
181               END DO
182            END DO
183         ENDIF ! ln_shelf_flx
184         DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
185            DO ji = 1, jpi
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(zwnd_i(ji,jj)*zwnd_i(ji,jj) + zwnd_j(ji,jj)*zwnd_j(ji,jj))
191                   cs = 0.63 + (0.066 * totwindspd) 
192                   utau(ji,jj) = cs * (rhoa/rau0) * zwnd_i(ji,jj) * totwindspd
193                   vtau(ji,jj) = cs * (rhoa/rau0) * zwnd_j(ji,jj) * totwindspd
194                ELSE
195                   utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
196                   vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
197                ENDIF
198                qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1) 
199                IF( ln_foam_flx .OR. ln_shelf_flx ) THEN 
200                   !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot)
201                   qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) 
202                   !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P
203                   emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1) 
204                ELSE
205                   qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 
206                   emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 
207                ENDIF
208               !!sfx (ji,jj) = sf(jp_sfx )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
209            END DO
210         END DO
211         !
212         IF( ln_shelf_flx ) THEN
213            ! calculate first the wind module, as it will be used later
214            DO jj = 2, jpjm1
215               DO ji = fs_2, fs_jpim1    ! vect. opt.
216                  ztx = zwnd_i(ji-1,jj  ) + zwnd_i(ji,jj)
217                  zty = zwnd_j(ji  ,jj-1) + zwnd_j(ji,jj)
218                  wndm(ji,jj) = 0.5 * SQRT( ztx * ztx + zty * zty )
219               END DO
220            END DO
221            CALL lbc_lnk_multi( 'sbcflx', wndm, 'T', 1. )
222         ENDIF
223         !                                                        ! add to qns the heat due to e-p
224         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
225         !
226         ! clem: without these lbc calls, it seems that the northfold is not ok (true in 3.6, not sure in 4.x)
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      !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe
246      IF( ln_foam_flx ) THEN
247         CALL lbc_lnk_multi( 'sbcflx', utau, 'U', -1._wp, vtau, 'V', -1._wp ) 
248      ENDIF
249      zcoef = 1. / ( zrhoa * zcdrag )
250      DO jj = 2, jpjm1
251         DO ji = fs_2, fs_jpim1   ! vect. opt.
252            ztx = ( utau(ji-1,jj  ) + utau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( umask(ji-1,jj  ,1), umask(ji,jj,1) ) )
253            zty = ( vtau(ji  ,jj-1) + vtau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( vmask(ji  ,jj-1,1), vmask(ji,jj,1) ) ) 
254            zmod = 0.5_wp * SQRT( ztx * ztx + zty * zty ) * tmask(ji,jj,1)
255            taum(ji,jj) = zmod
256            IF( .NOT. ln_shelf_flx ) THEN
257               wndm(ji,jj) = SQRT( zmod * zcoef )  !!clem: not used?
258            ENDIF
259         END DO
260      END DO
261      !
262      CALL lbc_lnk_multi( 'sbcflx', taum, 'T', 1._wp, wndm, 'T', 1._wp )
263      !
264      IF( ln_shelf_flx ) THEN
265         DEALLOCATE( zwnd_i, zwnd_j )
266      ENDIF
267      !
268   END SUBROUTINE sbc_flx
269
270   !!======================================================================
271END MODULE sbcflx
Note: See TracBrowser for help on using the repository browser.