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

Last change on this file since 15416 was 15416, checked in by hadjt, 12 months ago

SBC/sbcflx.F90

Fixing a bug in the shelf treatement of qsr and qns input.

File size: 13.8 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               IF( ln_shelf_flx ) THEN
194                  !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot)
195                  qns (ji,jj) = (sf(jp_qtot)%fnow(ji,jj,1)) * tmask(ji,jj,1)
196                  !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P
197                  emp (ji,jj) = (-1. * sf(jp_emp )%fnow(ji,jj,1)) * tmask(ji,jj,1)
198               ELSE
199                  qns (ji,jj) = (sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)) * tmask(ji,jj,1)
200                  emp (ji,jj) = (sf(jp_emp )%fnow(ji,jj,1)) * tmask(ji,jj,1)
201               ENDIF
202!               !utau(ji,jj) =   sf(jp_utau)%fnow(ji,jj,1)                              * umask(ji,jj,1)
203!               !vtau(ji,jj) =   sf(jp_vtau)%fnow(ji,jj,1)                              * vmask(ji,jj,1)
204!               !JT qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
205!               qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - qsr (ji,jj)              ) * tmask(ji,jj,1)
206!               emp (ji,jj) =   sf(jp_emp )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
207!               !!sfx (ji,jj) = sf(jp_sfx )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
208            END DO
209         END DO
210         !                                                        ! add to qns the heat due to e-p
211         !clem: I do not think it is needed
212         !!qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
213         !
214         ! clem: without these lbc calls, it seems that the northfold is not ok (true in 3.6, not sure in 4.x)
215
216         !JT Enabling as used in previous version using these climate forcings
217         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
218         !JT
219
220         CALL lbc_lnk_multi( 'sbcflx', utau, 'U', -1._wp, vtau, 'V', -1._wp, &
221            &                           qns, 'T',  1._wp, emp , 'T',  1._wp, qsr, 'T', 1._wp ) !! sfx, 'T', 1._wp  )
222         !
223         IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
224            WRITE(numout,*) 
225            WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK'
226            DO jf = 1, jpfld_local
227               IF( jf == jp_utau .OR. jf == jp_vtau )   zfact =     1.
228               IF( jf == jp_qtot .OR. jf == jp_qsr  )   zfact =     0.1
229               IF( jf == jp_emp                     )   zfact = 86400.
230               WRITE(numout,*) 
231               WRITE(numout,*) ' day: ', ndastp , TRIM(sf(jf)%clvar), ' * ', zfact
232            END DO
233         ENDIF
234         !
235      ENDIF
236      !                                                           ! module of wind stress and wind speed at T-point
237      ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
238      zcoef = 1. / ( zrhoa * zcdrag )
239      DO jj = 2, jpjm1
240         DO ji = fs_2, fs_jpim1   ! vect. opt.
241            ztx = ( utau(ji-1,jj  ) + utau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( umask(ji-1,jj  ,1), umask(ji,jj,1) ) )
242            zty = ( vtau(ji  ,jj-1) + vtau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( vmask(ji  ,jj-1,1), vmask(ji,jj,1) ) ) 
243            zmod = 0.5_wp * SQRT( ztx * ztx + zty * zty ) * tmask(ji,jj,1)
244            taum(ji,jj) = zmod
245            wndm(ji,jj) = SQRT( zmod * zcoef )  !!clem: not used?
246         END DO
247      END DO
248      !
249      CALL lbc_lnk_multi( 'sbcflx', taum, 'T', 1._wp, wndm, 'T', 1._wp )
250      !
251      !
252   END SUBROUTINE sbc_flx
253
254   !!======================================================================
255END MODULE sbcflx
Note: See TracBrowser for help on using the repository browser.