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 branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90 @ 7792

Last change on this file since 7792 was 7792, checked in by jcastill, 7 years ago

Changes as in branch branches/UKMO/r6232_INGV1_WAVE-coupling@7763
File size: 12.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   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 sbcwave         ! wave physics
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC sbc_flx       ! routine called by step.F90
30
31   INTEGER , PARAMETER ::   jpfld   = 5   ! maximum number of files to read
32   INTEGER , PARAMETER ::   jp_utau = 1   ! index of wind stress (i-component) file
33   INTEGER , PARAMETER ::   jp_vtau = 2   ! index of wind stress (j-component) file
34   INTEGER , PARAMETER ::   jp_qtot = 3   ! index of total (non solar+solar) heat file
35
36   ! These lines are here just to allow merging with another branch
37   INTEGER, PUBLIC            ::   nn_drag        ! type of formula to calculate wind stress from wind components
38   INTEGER, PUBLIC, PARAMETER ::   jp_ukmo  = 0        ! UKMO SHELF formulation
39   INTEGER, PUBLIC, PARAMETER ::   jp_std   = 1        ! standard formulation with forced or coupled drag coefficient
40   INTEGER, PUBLIC, PARAMETER ::   jp_const = 2        ! standard formulation with constant drag coefficient
41
42   INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file
43   INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file
44   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
45   LOGICAL , PUBLIC    ::   ln_shelf_flx = .FALSE. ! UKMO SHELF specific flux flag
46   INTEGER             ::   jpfld_local   ! maximum number of files to read (locally modified depending on ln_shelf_flx)
47
48   !! * Substitutions
49#  include "domzgr_substitute.h90"
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/OPA 3.3 , NEMO-consortium (2010)
53   !! $Id$
54   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
55   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE sbc_flx( kt )
59      !!---------------------------------------------------------------------
60      !!                    ***  ROUTINE sbc_flx  ***
61      !!                   
62      !! ** Purpose :   provide at each time step the surface ocean fluxes
63      !!                (momentum, heat, freshwater and runoff)
64      !!
65      !! ** Method  : - READ each fluxes in NetCDF files:
66      !!                   i-component of the stress              utau  (N/m2)
67      !!                   j-component of the stress              vtau  (N/m2)
68      !!                   net downward heat flux                 qtot  (watt/m2)
69      !!                   net downward radiative flux            qsr   (watt/m2)
70      !!                   net upward freshwater (evapo - precip) emp   (kg/m2/s)
71      !!
72      !!      CAUTION :  - never mask the surface stress fields
73      !!                 - the stress is assumed to be in the (i,j) mesh referential
74      !!
75      !! ** Action  :   update at each time-step
76      !!              - utau, vtau  i- and j-component of the wind stress
77      !!              - taum        wind stress module at T-point
78      !!              - wndm        10m wind module at T-point
79      !!              - qns         non solar heat flux including heat flux due to emp
80      !!              - qsr         solar heat flux
81      !!              - emp         upward mass flux (evap. - precip.)
82      !!              - sfx         salt flux; set to zero at nit000 but possibly non-zero
83      !!                            if ice is present (computed in limsbc(_2).F90)
84      !!----------------------------------------------------------------------
85      INTEGER, INTENT(in) ::   kt   ! ocean time step
86      !!
87      INTEGER  ::   ji, jj, jf            ! dummy indices
88      INTEGER  ::   ierror                ! return error code
89      INTEGER  ::   ios                   ! Local integer output status for namelist read
90      REAL(wp) ::   zfact                 ! temporary scalar
91      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
92      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
93      REAL(wp) ::   totwind               ! UKMO SHELF: Module of wind speed
94      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
95      !!
96      CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files
97      NAMELIST/namsbc_flx/ ln_shelf_flx                           ! Put here to allow merging with another UKMO branch
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  ! informations about the fields to be read
100      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp
101      !!---------------------------------------------------------------------
102      !
103      IF( kt == nit000 ) THEN                ! First call kt=nit000 
104         ! set file information
105         REWIND( numnam_ref )              ! Namelist namsbc_flx in reference namelist : Files for fluxes
106         READ  ( numnam_ref, namsbc_flx, IOSTAT = ios, ERR = 901)
107901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_flx in reference namelist', lwp )
108
109         REWIND( numnam_cfg )              ! Namelist namsbc_flx in configuration namelist : Files for fluxes
110         READ  ( numnam_cfg, namsbc_flx, IOSTAT = ios, ERR = 902 )
111902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_flx in configuration namelist', lwp )
112         IF(lwm) WRITE ( numond, namsbc_flx ) 
113         !
114         !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing?
115         IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 )   &
116            &   CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' ) 
117         !
118         !                                         ! store namelist information in an array
119         slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau
120         slf_i(jp_qtot) = sn_qtot   ;   slf_i(jp_qsr ) = sn_qsr 
121         slf_i(jp_emp ) = sn_emp
122         !
123         ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
124         IF( ierror > 0 ) THEN   
125            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN 
126         ENDIF
127         DO ji= 1, jpfld
128            ALLOCATE( sf(ji)%fnow(jpi,jpj,1) )
129            IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) )
130         END DO
131         !                                         ! fill sf with slf_i and control print
132         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
133         !
134         sfx(:,:) = 0.0_wp                         ! salt flux due to freezing/melting (non-zero only if ice is present; set in limsbc(_2).F90)
135         !
136      ENDIF
137
138      CALL fld_read( kt, nn_fsbc, sf )                            ! input fields provided at the current time-step
139     
140      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency
141
142         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )   ! modify now Qsr to include the diurnal cycle
143         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)
144         ENDIF
145!CDIR COLLAPSE
146         DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
147            DO ji = 1, jpi
148               utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
149               vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
150               qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)
151               emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1)
152            END DO
153         END DO
154         !                                                        ! add modification due to drag coefficient read from wave forcing
155         !                                                        ! this code is inefficient but put here to allow merging with another UKMO branch
156         IF( ln_shelf_flx ) THEN
157            IF( ln_cdgw .AND. nn_drag == jp_std ) THEN
158               IF( cpl_wdrag ) THEN 
159                  ! reset utau and vtau to the wind components: the momentum will
160                  ! be calculated from the coupled value of the drag coefficient
161                  DO jj = 1, jpj
162                     DO ji = 1, jpi
163                        utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
164                        vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
165                     END DO
166                  END DO
167               ELSE
168                  DO jj = 1, jpj
169                     DO ji = 1, jpi
170                        totwind = sqrt((sf(jp_utau)%fnow(ji,jj,1))**2.0 + (sf(jp_vtau)%fnow(ji,jj,1))**2.0)
171                        utau(ji,jj) = zrhoa * cdn_wave(ji,jj) * sf(jp_utau)%fnow(ji,jj,1) * totwind
172                        vtau(ji,jj) = zrhoa * cdn_wave(ji,jj) * sf(jp_vtau)%fnow(ji,jj,1) * totwind
173                     END DO
174                  END DO
175               ENDIF
176            ELSE IF( nn_drag == jp_const ) THEN
177               DO jj = 1, jpj
178                  DO ji = 1, jpi
179                     totwind = sqrt((sf(jp_utau)%fnow(ji,jj,1))**2.0 + (sf(jp_vtau)%fnow(ji,jj,1))**2.0)
180                     utau(ji,jj) = zrhoa * zcdrag * sf(jp_utau)%fnow(ji,jj,1) * totwind
181                     vtau(ji,jj) = zrhoa * zcdrag * sf(jp_vtau)%fnow(ji,jj,1) * totwind
182                  END DO
183               END DO
184            ENDIF
185         ENDIF
186         !                                                        ! add to qns the heat due to e-p
187         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
188         !
189         !                                                        ! module of wind stress and wind speed at T-point
190         zcoef = 1. / ( zrhoa * zcdrag )
191!CDIR NOVERRCHK
192         DO jj = 2, jpjm1
193!CDIR NOVERRCHK
194            DO ji = fs_2, fs_jpim1   ! vect. opt.
195               ztx = utau(ji-1,jj  ) + utau(ji,jj) 
196               zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
197               zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
198               taum(ji,jj) = zmod
199               IF( ln_shelf_flx ) THEN
200                  ztx = sf(jp_utau)%fnow(ji-1,jj  ) + sf(jp_utau)%fnow(ji,jj)
201                  zty = sf(jp_vtau)%fnow(ji  ,jj-1) + sf(jp_vtau)%fnow(ji,jj)
202                  wndm(ji,jj) = 0.5 * SQRT( ztx * ztx + zty * zty )
203               ELSE
204               wndm(ji,jj) = SQRT( zmod * zcoef )
205               ENDIF
206            END DO
207         END DO
208         taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1)
209         CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
210
211         IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
212            WRITE(numout,*) 
213            WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK'
214            DO jf = 1, jpfld
215               IF( jf == jp_utau .OR. jf == jp_vtau )   zfact =     1.
216               IF( jf == jp_qtot .OR. jf == jp_qsr  )   zfact =     0.1
217               IF( jf == jp_emp                     )   zfact = 86400.
218               WRITE(numout,*) 
219               WRITE(numout,*) ' day: ', ndastp , TRIM(sf(jf)%clvar), ' * ', zfact
220               CALL prihre( sf(jf)%fnow, jpi, jpj, 1, jpi, 20, 1, jpj, 10, zfact, numout )
221            END DO
222            CALL FLUSH(numout)
223         ENDIF
224         !
225      ENDIF
226      !
227   END SUBROUTINE sbc_flx
228
229   !!======================================================================
230END MODULE sbcflx
Note: See TracBrowser for help on using the repository browser.