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/releases/r4.0/r4.0-HEAD/src/OCE/SBC – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcflx.F90 @ 15813

Last change on this file since 15813 was 15813, checked in by clem, 8 months ago

SI3: change heat budget to avoid supercooling. The rest: cosmetics
  • Property svn:keywords set to Id
File size: 12.2 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 trc_oce         ! share SMS/Ocean variables
18   USE sbcdcy          ! surface boundary condition: diurnal cycle on qsr
19   USE phycst          ! physical constants
20   !
21   USE fldread         ! read input fields
22   USE iom             ! IOM library
23   USE in_out_manager  ! I/O manager
24   USE lib_mpp         ! distribued memory computing library
25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC sbc_flx       ! routine called by step.F90
31
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   INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file
36   INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file
37!!$   INTEGER , PARAMETER ::   jp_sfx  = 6   ! index of salt flux
38!!$   INTEGER , PARAMETER ::   jp_sithic = 7    ! index of sea ice thickness
39!!$   INTEGER , PARAMETER ::   jp_siconc = 8    ! index of sea ice fraction
40!!$   INTEGER , PARAMETER ::   jp_hfrnf  = 9    ! index of rnf heat flux
41!!$   INTEGER , PARAMETER ::   jp_hfisf  = 10   ! index of iceshelf heat flux
42!!$   INTEGER , PARAMETER ::   jp_fwisf  = 11   ! index of iceshelf freshwater flux
43!!$   INTEGER , PARAMETER ::   jp_fwrnf  = 12   ! index of runoff freshwater flux
44   INTEGER , PARAMETER ::   jpfld   = 5   ! maximum number of files to read
45
46   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
47
48   !! * Substitutions
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
52   !! $Id$
53   !! Software governed by the CeCILL license (see ./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE sbc_flx( kt )
58      !!---------------------------------------------------------------------
59      !!                    ***  ROUTINE sbc_flx  ***
60      !!                   
61      !! ** Purpose :   provide at each time step the surface ocean fluxes
62      !!                (momentum, heat, freshwater and runoff)
63      !!
64      !! ** Method  : - READ each fluxes in NetCDF files:
65      !!                   i-component of the stress              utau  (N/m2)
66      !!                   j-component of the stress              vtau  (N/m2)
67      !!                   net downward heat flux                 qtot  (watt/m2)
68      !!                   net downward radiative flux            qsr   (watt/m2)
69      !!                   net upward freshwater (evapo - precip) emp   (kg/m2/s)
70      !!                   salt flux                              sfx   (pss*dh*rho/dt => g/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 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      !!
94      CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files
95      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                    ! array of namelist information structures
96      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp  !!, & ! informations about the fields to be read
97!!$         &             sn_sfx, sn_sithic, sn_siconc, sn_hfisf, sn_hfrnf, sn_fwisf, sn_fwrnf
98      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp !!, &
99!!$         &                         sn_sfx, sn_sithic, sn_siconc, sn_hfisf, sn_hfrnf, sn_fwisf, sn_fwrnf
100      !!---------------------------------------------------------------------
101      !
102      IF( kt == nit000 ) THEN                ! First call kt=nit000 
103         ! set file information
104         REWIND( numnam_ref )              ! Namelist namsbc_flx in reference namelist : Files for fluxes
105         READ  ( numnam_ref, namsbc_flx, IOSTAT = ios, ERR = 901)
106901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_flx in reference namelist' )
107
108         REWIND( numnam_cfg )              ! Namelist namsbc_flx in configuration namelist : Files for fluxes
109         READ  ( numnam_cfg, namsbc_flx, IOSTAT = ios, ERR = 902 )
110902      IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_flx in configuration namelist' )
111         IF(lwm) WRITE ( numond, namsbc_flx ) 
112         !
113         !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing?
114         IF( ln_dm2dc .AND. sn_qsr%freqh /= 24. )   &
115            &   CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' ) 
116         !
117         !                                         ! store namelist information in an array
118         slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau
119         slf_i(jp_qtot) = sn_qtot   ;   slf_i(jp_qsr ) = sn_qsr 
120         slf_i(jp_emp ) = sn_emp
121!!$         slf_i(jp_sfx  ) = sn_sfx
122!!$         slf_i(jp_sithic) = sn_sithic
123!!$         slf_i(jp_siconc) = sn_siconc
124!!$         slf_i(jp_hfisf) = sn_hfisf    ;   slf_i(jp_hfrnf) = sn_hfrnf
125!!$         slf_i(jp_fwisf) = sn_fwisf    ;   slf_i(jp_fwrnf) = sn_fwrnf
126         !
127         ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
128         IF( ierror > 0 ) THEN   
129            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN 
130         ENDIF
131         DO ji= 1, jpfld
132            ALLOCATE( sf(ji)%fnow(jpi,jpj,1) )
133            IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) )
134         END DO
135         !                                         ! fill sf with slf_i and control print
136         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
137         !
138      ENDIF
139
140      CALL fld_read( kt, nn_fsbc, sf )                            ! input fields provided at the current time-step
141     
142      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency
143
144         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)  ! modify now Qsr to include the diurnal cycle
145         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
146         ENDIF
147#if defined key_top
148      IF( ln_trcdc2dm )  THEN      !  diurnal cycle in TOP
149         IF( ln_dm2dc )  THEN  ;  qsr_mean(:,:) = sf(jp_qsr)%fnow(:,:,1)  * tmask(:,:,1)
150         ELSE                  ;  ncpl_qsr_freq = sf(jp_qsr)%freqh * 3600 !  qsr_mean will be computed in TOP
151         ENDIF
152      ENDIF
153#endif         
154         DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
155            DO ji = 1, jpi
156               utau(ji,jj) =   sf(jp_utau)%fnow(ji,jj,1)                              * umask(ji,jj,1)
157               vtau(ji,jj) =   sf(jp_vtau)%fnow(ji,jj,1)                              * vmask(ji,jj,1)
158               qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
159               emp (ji,jj) =   sf(jp_emp )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
160!!$               sfx (ji,jj) = sf(jp_sfx   )%fnow(ji,jj,1)                             * tmask(ji,jj,1)
161!!$               !! => if the following is used, then one needs to change tke routine + allocate hm_i in sbc_oce
162!!$               hm_i(ji,jj) = sf(jp_sithic)%fnow(ji,jj,1)                             * tmask(ji,jj,1)
163!!$               fr_i(ji,jj) = sf(jp_siconc)%fnow(ji,jj,1)                             * tmask(ji,jj,1)
164!!$               !! => if the following is used, then one needs to change rnf and isf routines + allocate the arrays
165!!$               hfisf(ji,jj) = sf(jp_hfisf)%fnow(ji,jj,1)                             * ssmask(ji,jj)
166!!$               fwisf(ji,jj) = sf(jp_fwisf)%fnow(ji,jj,1)                             * ssmask(ji,jj)
167!!$               hfrnf(ji,jj) = sf(jp_hfrnf)%fnow(ji,jj,1)                             * tmask(ji,jj,1)
168!!$               fwrnf(ji,jj) = sf(jp_fwrnf)%fnow(ji,jj,1)                             * tmask(ji,jj,1)
169            END DO
170         END DO
171         !
172         ! clem: without these lbc calls, it seems that the northfold is not ok (true in 3.6, not sure in 4.x)
173         CALL lbc_lnk_multi( 'sbcflx', utau, 'U', -1._wp, vtau, 'V', -1._wp, &
174            &                           qns, 'T',  1._wp, emp , 'T',  1._wp, qsr, 'T', 1._wp ) !! &
175!!$            &                           sfx, 'T', 1._wp, hm_i,'T', 1._wp, fr_i,'T', 1._wp, &
176!!$            &                          hfisf, 'T', 1._wp, fwisf, 'T', 1._wp, hfrnf, 'T', 1._wp, fwrnf, 'T', 1._wp )
177         !
178         IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
179            WRITE(numout,*) 
180            WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK'
181            DO jf = 1, jpfld
182               IF( jf == jp_utau .OR. jf == jp_vtau )   zfact =     1.
183               IF( jf == jp_qtot .OR. jf == jp_qsr  )   zfact =     0.1
184               IF( jf == jp_emp                     )   zfact = 86400.
185               WRITE(numout,*) 
186               WRITE(numout,*) ' day: ', ndastp , TRIM(sf(jf)%clvar), ' * ', zfact
187            END DO
188         ENDIF
189         !
190      ENDIF
191      !                                                           ! module of wind stress and wind speed at T-point
192      ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
193      zcoef = 1. / ( zrhoa * zcdrag )
194      DO jj = 2, jpjm1
195         DO ji = fs_2, fs_jpim1   ! vect. opt.
196            ztx = ( utau(ji-1,jj  ) + utau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( umask(ji-1,jj  ,1), umask(ji,jj,1) ) )
197            zty = ( vtau(ji  ,jj-1) + vtau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( vmask(ji  ,jj-1,1), vmask(ji,jj,1) ) ) 
198            zmod = SQRT( ztx * ztx + zty * zty ) * tmask(ji,jj,1)
199            taum(ji,jj) = zmod
200            wndm(ji,jj) = SQRT( zmod * zcoef )  !!clem: not used?
201         END DO
202      END DO
203      !
204      CALL lbc_lnk_multi( 'sbcflx', taum, 'T', 1._wp, wndm, 'T', 1._wp )
205      !
206      !
207   END SUBROUTINE sbc_flx
208
209   !!======================================================================
210END MODULE sbcflx
Note: See TracBrowser for help on using the repository browser.