source: branches/UKMO/AMM15_v3_6_STABLE_package_UKEP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90 @ 8757

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

Changes to the branch to make it merge with branch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/r6232_HZG_WAVE-coupling
WARNING: This branch will probably not compile on its own - changes in the original branch will have to be copied here

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