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

Last change on this file since 11277 was 11277, checked in by kingr, 20 months ago

Merged Juan's changes for running AMM15 woth wave coupling.
Corrected minor logic error to allow AMM7-uncoupled to reproduce earlier results.
Few line spacing changes to allow merging with OBS br and trunk rvn 5518.

File size: 17.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   USE fldread         ! read input fields
20   USE iom             ! IOM library
21   USE in_out_manager  ! I/O manager
22   USE sbcwave         ! wave physics
23   USE lib_mpp         ! distribued memory computing library
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE wrk_nemo        ! work arrays
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC sbc_flx       ! routine called by step.F90
31
32   INTEGER , PARAMETER ::   jpfld   = 6   ! maximum number of files to read
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
38   INTEGER , PARAMETER ::   jp_press = 6  ! index of pressure for UKMO shelf fluxes
39   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
40   LOGICAL , PUBLIC    ::   ln_shelf_flx = .FALSE. ! UKMO SHELF specific flux flag
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)
43   INTEGER             ::   jpfld_local   ! maximum number of files to read (locally modified depending on ln_shelf_flx)
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.3 , NEMO-consortium (2010)
50   !! $Id$
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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
70      !!                 - the stress is assumed to be in the (i,j) mesh referential
71      !!
72      !! ** Action  :   update at each time-step
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
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)
81      !!----------------------------------------------------------------------
82      INTEGER, INTENT(in) ::   kt   ! ocean time step
83      !!
84      INTEGER  ::   ji, jj, jf            ! dummy indices
85      INTEGER  ::   ierror                ! return error code
86      INTEGER  ::   ios                   ! Local integer output status for namelist read
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
90      REAL(wp) ::   totwind               ! UKMO SHELF: Module of wind speed
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), DIMENSION(:,:), POINTER ::   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', 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 )
117         IF(lwm) WRITE ( numond, namsbc_flx ) 
118         !
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
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
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
136         !
137            ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
138            IF( ln_shelf_flx ) slf_i(jp_press) = sn_press
139   
140            ! define local jpfld depending on shelf_flx logical
141            IF( ln_shelf_flx ) THEN
142               jpfld_local = jpfld
143            ELSE
144               jpfld_local = jpfld-1
145            ENDIF
146            !
147         IF( ierror > 0 ) THEN   
148            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN 
149         ENDIF
150         DO ji= 1, jpfld_local
151            ALLOCATE( sf(ji)%fnow(jpi,jpj,1) )
152            IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) )
153         END DO
154         !                                         ! fill sf with slf_i and control print
155         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
156         !
157         sfx(:,:) = 0.0_wp                         ! salt flux due to freezing/melting (non-zero only if ice is present; set in limsbc(_2).F90)
158         !
159      ENDIF
160
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
164
165         !!UKMO SHELF wind speed relative to surface currents - put here to allow merging with coupling branch
166         IF( ln_shelf_flx ) THEN
167            CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j )
168
169            IF( ln_rel_wind ) THEN
170               DO jj = 1, jpj
171                  DO ji = 1, jpi
172                     zwnd_i(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) - rn_wfac * ssu_m(ji,jj)
173                     zwnd_j(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) - rn_wfac * ssv_m(ji,jj)
174                  END DO
175               END DO
176            ELSE
177               zwnd_i(:,:) = sf(jp_utau)%fnow(:,:,1)
178               zwnd_j(:,:) = sf(jp_vtau)%fnow(:,:,1)
179            ENDIF
180         ENDIF
181
182         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )   ! modify now Qsr to include the diurnal cycle
183         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)
184         ENDIF
185!CDIR COLLAPSE
186            !!UKMO SHELF effect of atmospheric pressure on SSH
187            ! If using ln_apr_dyn, this is done there so don't repeat here.
188            IF( ln_shelf_flx .AND. .NOT. ln_apr_dyn) THEN
189               DO jj = 1, jpjm1
190                  DO ji = 1, jpim1
191                     apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj)
192                     apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj)
193                  END DO
194               END DO
195            ENDIF ! ln_shelf_flx
196
197         DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
198            DO ji = 1, jpi
199                   IF( ln_shelf_flx ) THEN
200                      !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing
201                      pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1)
202                      !! UKMO SHELF flux files contain wind speed not wind stress
203                      totwindspd = sqrt(zwnd_i(ji,jj)*zwnd_i(ji,jj) + zwnd_j(ji,jj)*zwnd_j(ji,jj))
204                      cs = 0.63 + (0.066 * totwindspd)
205                      utau(ji,jj) = cs * (rhoa/rau0) * zwnd_i(ji,jj) * totwindspd
206                      vtau(ji,jj) = cs * (rhoa/rau0) * zwnd_j(ji,jj) * totwindspd
207                   ELSE
208                      utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
209                      vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
210                   ENDIF
211                   qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1)
212                   IF( ln_foam_flx .OR. ln_shelf_flx ) THEN
213                      !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot)
214                      qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1)
215                      !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P
216                      emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1)
217                   ELSE
218                      qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)
219                      emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1)
220                   ENDIF
221            END DO
222         END DO
223         !                                                        ! add modification due to drag coefficient read from wave forcing
224         !                                                        ! this code is inefficient but put here to allow merging with another UKMO branch
225         IF( ln_shelf_flx ) THEN 
226            ! calculate first the wind module, as it will be used later
227            DO jj = 2, jpjm1 
228               DO ji = fs_2, fs_jpim1    ! vect. opt.
229                  ztx = zwnd_i(ji-1,jj  ) + zwnd_i(ji,jj) 
230                  zty = zwnd_j(ji  ,jj-1) + zwnd_j(ji,jj) 
231                  wndm(ji,jj) = 0.5 * SQRT( ztx * ztx + zty * zty ) 
232               END DO
233            END DO
234            CALL lbc_lnk( wndm(:,:), 'T', 1. ) 
235       
236            IF( ln_cdgw .AND. nn_drag == jp_std ) THEN
237               IF( cpl_wdrag ) THEN 
238                  ! reset utau and vtau to the wind components: the momentum will
239                  ! be calculated from the coupled value of the drag coefficient
240                  DO jj = 1, jpj 
241                     DO ji = 1, jpi 
242                        utau(ji,jj) = zwnd_i(ji,jj) 
243                        vtau(ji,jj) = zwnd_j(ji,jj) 
244                     END DO
245                  END DO
246               ELSE
247                  DO jj = 1, jpjm1 
248                     DO ji = 1, jpim1 
249                        utau(ji,jj) = zrhoa * 0.5 * ( cdn_wave(ji,jj) + cdn_wave(ji+1,jj) ) * zwnd_i(ji,jj) * & 
250                                                                        0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) ) 
251                        vtau(ji,jj) = zrhoa * 0.5 * ( cdn_wave(ji,jj) + cdn_wave(ji,jj+1) ) * zwnd_j(ji,jj) * & 
252                                                                        0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) ) 
253                     END DO
254                  END DO
255                  CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. ) 
256               ENDIF
257            ELSE IF( nn_drag == jp_const ) THEN
258               DO jj = 1, jpjm1 
259                  DO ji = 1, jpim1 
260                     utau(ji,jj) = zrhoa * zcdrag * zwnd_i(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) ) 
261                     vtau(ji,jj) = zrhoa * zcdrag * zwnd_j(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) ) 
262                  END DO
263               END DO
264               CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. ) 
265            ENDIF
266         ENDIF
267         !                                                        ! add to qns the heat due to e-p
268         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
269         !
270   
271            !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe
272            IF( ln_foam_flx ) THEN
273               CALL lbc_lnk( utau(:,:), 'U', -1. )
274               CALL lbc_lnk( vtau(:,:), 'V', -1. )
275            ENDIF
276   
277         !                                                        ! module of wind stress and wind speed at T-point
278         zcoef = 1. / ( zrhoa * zcdrag )
279!CDIR NOVERRCHK
280         DO jj = 2, jpjm1
281!CDIR NOVERRCHK
282            DO ji = fs_2, fs_jpim1   ! vect. opt.
283               ztx = utau(ji-1,jj  ) + utau(ji,jj) 
284               zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
285               zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
286               taum(ji,jj) = zmod
287               IF ( .NOT. (ln_shelf_flx .AND. ln_cpl)) THEN
288               wndm(ji,jj) = SQRT( zmod * zcoef )
289               ENDIF
290            END DO
291         END DO
292         taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1)
293         CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
294
295         IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
296            WRITE(numout,*) 
297            WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK'
298            DO jf = 1, jpfld_local
299               IF( jf == jp_utau .OR. jf == jp_vtau )   zfact =     1.
300               IF( jf == jp_qtot .OR. jf == jp_qsr  )   zfact =     0.1
301               IF( jf == jp_emp                     )   zfact = 86400.
302               WRITE(numout,*) 
303               WRITE(numout,*) ' day: ', ndastp , TRIM(sf(jf)%clvar), ' * ', zfact
304               CALL prihre( sf(jf)%fnow, jpi, jpj, 1, jpi, 20, 1, jpj, 10, zfact, numout )
305            END DO
306            CALL FLUSH(numout)
307         ENDIF
308         !
309         IF( ln_shelf_flx ) THEN
310            CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j )
311         ENDIF
312         !
313      ENDIF
314      !
315   END SUBROUTINE sbc_flx
316
317   !!======================================================================
318END MODULE sbcflx
Note: See TracBrowser for help on using the repository browser.