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.
sbcana.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta

  • Property svn:keywords set to Id
File size: 15.6 KB
Line 
1MODULE sbcana
2   !!======================================================================
3   !!                       ***  MODULE  sbcana  ***
4   !! Ocean forcing:  analytical momentum, heat and freshwater forcings
5   !!=====================================================================
6   !! History :  3.0   ! 2006-06  (G. Madec)  Original code
7   !!            3.2   ! 2009-07  (G. Madec)  Style only
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   sbc_ana  : set an analytical ocean forcing
12   !!   sbc_gyre : set the GYRE configuration analytical forcing
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 phycst          ! physical constants
18   USE in_out_manager  ! I/O manager
19   USE lib_mpp         ! distribued memory computing library
20   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   sbc_ana    ! routine called in sbcmod module
26   PUBLIC   sbc_gyre   ! routine called in sbcmod module
27
28   !                               !!* Namelist namsbc_ana *
29   INTEGER  ::   nn_tau000 = 1      ! nb of time-step during which the surface stress
30   !                                ! increase from 0 to its nominal value
31   REAL(wp) ::   rn_utau0  = 0.e0   ! constant wind stress value in i-direction
32   REAL(wp) ::   rn_vtau0  = 0.e0   ! constant wind stress value in j-direction
33   REAL(wp) ::   rn_qns0   = 0.e0   ! non solar heat flux
34   REAL(wp) ::   rn_qsr0   = 0.e0   !     solar heat flux
35   REAL(wp) ::   rn_emp0   = 0.e0   ! net freshwater flux
36   
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
42   !! $Id$
43   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45
46CONTAINS
47
48   SUBROUTINE sbc_ana( kt )
49      !!---------------------------------------------------------------------
50      !!                    ***  ROUTINE sbc_ana ***
51      !!             
52      !! ** Purpose :   provide at each time-step the ocean surface boundary
53      !!              condition, i.e. the momentum, heat and freshwater fluxes.
54      !!
55      !! ** Method  :   Constant and uniform surface forcing specified from
56      !!              namsbc_ana namelist parameters. All the fluxes are time
57      !!              independant except the stresses which increase from zero
58      !!              during the first nn_tau000 time-step
59      !!                CAUTION : never mask the surface stress field !
60      !!
61      !! ** Action  : - set the ocean surface boundary condition, i.e. 
62      !!                   utau, vtau, taum, wndm, qns, qsr, emp, emps
63      !!----------------------------------------------------------------------
64      INTEGER, INTENT(in) ::   kt       ! ocean time step
65      !!
66      REAL(wp)            ::   zfacto                ! local scalar
67      REAL(wp)            ::   zrhoa  = 1.22         ! Air density kg/m3
68      REAL(wp)            ::   zcdrag = 1.5e-3       ! drag coefficient
69      REAL(wp)            ::   ztx, zty, zmod, zcoef ! temporary variables
70      !!
71      NAMELIST/namsbc_ana/ nn_tau000, rn_utau0, rn_vtau0, rn_qns0, rn_qsr0, rn_emp0
72      !!---------------------------------------------------------------------
73      !
74      IF( kt == nit000 ) THEN
75         !
76         REWIND ( numnam )                   ! Read Namelist namsbc : surface fluxes
77         READ   ( numnam, namsbc_ana )
78         !
79         IF(lwp) WRITE(numout,*)' '
80         IF(lwp) WRITE(numout,*)' sbc_ana : Constant surface fluxes read in namsbc_ana namelist'
81         IF(lwp) WRITE(numout,*)' ~~~~~~~ '
82         IF(lwp) WRITE(numout,*)'              spin up of the stress  nn_tau000 = ', nn_tau000, ' time-steps'
83         IF(lwp) WRITE(numout,*)'              constant i-stress      rn_utau0  = ', rn_utau0 , ' N/m2'
84         IF(lwp) WRITE(numout,*)'              constant j-stress      rn_vtau0  = ', rn_vtau0 , ' N/m2'
85         IF(lwp) WRITE(numout,*)'              non solar heat flux    rn_qns0   = ', rn_qns0  , ' W/m2'
86         IF(lwp) WRITE(numout,*)'              solar heat flux        rn_qsr0   = ', rn_qsr0  , ' W/m2'
87         IF(lwp) WRITE(numout,*)'              net heat flux          rn_emp0   = ', rn_emp0  , ' Kg/m2/s'
88         !
89         nn_tau000 = MAX( nn_tau000, 1 )   ! must be >= 1
90         !
91      ENDIF
92
93      qns   (:,:) = rn_qns0
94      qsr   (:,:) = rn_qsr0
95      emp   (:,:) = rn_emp0
96      emps  (:,:) = rn_emp0
97   
98      ! Increase the surface stress to its nominal value during the first nn_tau000 time-steps
99      IF( kt <= nn_tau000 ) THEN
100         zfacto = 0.5 * (  1. - COS( rpi * FLOAT( kt ) / FLOAT( nn_tau000 ) )  )
101         zcoef = 1. / ( zrhoa * zcdrag ) 
102         ztx = zfacto * rn_utau0
103         zty = zfacto * rn_vtau0
104         zmod = SQRT( ztx * ztx + zty * zty )
105         utau(:,:) = ztx
106         vtau(:,:) = zty
107         taum(:,:) = zmod
108         zmod = SQRT( zmod * zcoef )
109         wndm(:,:) = zmod
110      ENDIF
111      !
112   END SUBROUTINE sbc_ana
113
114
115   SUBROUTINE sbc_gyre( kt )
116      !!---------------------------------------------------------------------
117      !!                    ***  ROUTINE sbc_ana ***
118      !!             
119      !! ** Purpose :   provide at each time-step the GYRE surface boundary
120      !!              condition, i.e. the momentum, heat and freshwater fluxes.
121      !!
122      !! ** Method  :   analytical seasonal cycle for GYRE configuration.
123      !!                CAUTION : never mask the surface stress field !
124      !!
125      !! ** Action  : - set the ocean surface boundary condition, i.e.   
126      !!                   utau, vtau, taum, wndm, qns, qsr, emp, emps
127      !!
128      !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000.
129      !!----------------------------------------------------------------------
130      INTEGER, INTENT(in) ::   kt          ! ocean time step
131      !!
132      INTEGER  ::   ji, jj                 ! dummy loop indices
133      INTEGER  ::   zyear0                 ! initial year
134      INTEGER  ::   zmonth0                ! initial month
135      INTEGER  ::   zday0                  ! initial day
136      INTEGER  ::   zday_year0             ! initial day since january 1st
137      REAL(wp) ::   ztau     , ztau_sais   ! wind intensity and of the seasonal cycle
138      REAL(wp) ::   ztime                  ! time in hour
139      REAL(wp) ::   ztimemax , ztimemin    ! 21th June, and 21th decem. if date0 = 1st january
140      REAL(wp) ::   ztimemax1, ztimemin1   ! 21th June, and 21th decem. if date0 = 1st january
141      REAL(wp) ::   ztimemax2, ztimemin2   ! 21th June, and 21th decem. if date0 = 1st january
142      REAL(wp) ::   ztaun                  ! intensity
143      REAL(wp) ::   zemp_s, zemp_n, zemp_sais, ztstar
144      REAL(wp) ::   zcos_sais1, zcos_sais2, ztrp, zconv, t_star
145      REAL(wp) ::   zsumemp, zsurf
146      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
147      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
148      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
149      REAL(wp) ::   zyydd                 ! number of days in one year
150      !!---------------------------------------------------------------------
151      zyydd = REAL(nyear_len(1),wp)
152
153      ! ---------------------------- !
154      !  heat and freshwater fluxes  !
155      ! ---------------------------- !
156      !same temperature, E-P as in HAZELEGER 2000
157
158      zyear0     =   ndate0 / 10000                             ! initial year
159      zmonth0    = ( ndate0 - zyear0 * 10000 ) / 100            ! initial month
160      zday0      =   ndate0 - zyear0 * 10000 - zmonth0 * 100    ! initial day betwen 1 and 30
161      zday_year0 = ( zmonth0 - 1 ) * 30.+zday0                  ! initial day betwen 1 and 360
162
163      ! current day (in hours) since january the 1st of the current year
164      ztime = REAL( kt ) * rdt / (rmmss * rhhmm)   &       !  total incrementation (in hours)
165         &      - (nyear  - 1) * rjjhh * zyydd             !  minus years since beginning of experiment (in hours)
166
167      ztimemax1 = ((5.*30.)+21.)* 24.                      ! 21th june     at 24h in hours
168      ztimemin1 = ztimemax1 + rjjhh * zyydd / 2            ! 21th december        in hours
169      ztimemax2 = ((6.*30.)+21.)* 24.                      ! 21th july     at 24h in hours
170      ztimemin2 = ztimemax2 - rjjhh * zyydd / 2            ! 21th january         in hours
171      !                                                    ! NB: rjjhh * zyydd / 4 = one seasonal cycle in hours
172
173      ! amplitudes
174      zemp_S    = 0.7       ! intensity of COS in the South
175      zemp_N    = 0.8       ! intensity of COS in the North
176      zemp_sais = 0.1
177      zTstar    = 28.3      ! intemsity from 28.3 a -5 deg
178
179      ! 1/2 period between 21th June and 21th December and between 21th July and 21th January
180      zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi ) 
181      zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi )
182
183      ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K)
184      zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s
185      DO jj = 1, jpj
186         DO ji = 1, jpi
187            ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
188            ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
189            ! 64.5 in summer, 42.5 in winter
190            t_star = zTstar * ( 1 + 1. / 50. * zcos_sais2 )                &
191               &                    * COS( rpi * (gphit(ji,jj) - 5.)               &
192               &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) )
193            ! 23.5 deg : tropics
194            qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
195            qns (ji,jj) = ztrp * ( tb(ji,jj,1) - t_star ) - qsr(ji,jj)
196            IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg
197               emp  (ji,jj) =   zemp_S * zconv   &
198                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  &
199                  &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
200            ELSE
201               emp (ji,jj) =  - zemp_N * zconv   &
202                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  &
203                  &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
204            ENDIF
205         END DO
206      END DO
207      emps(:,:) = emp(:,:)
208
209      ! Compute the emp flux such as its integration on the whole domain at each time is zero
210      IF( nbench /= 1 ) THEN
211         zsumemp = 0.e0   ;   zsurf = 0.e0
212         DO jj = 1, jpj
213            DO ji = 1, jpi
214               zsumemp = zsumemp + emp(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
215               zsurf   = zsurf   +              tmask(ji,jj,1) * tmask_i(ji,jj)
216            END DO
217         END DO
218
219         IF( lk_mpp )   CALL mpp_sum( zsumemp  )       ! sum over the global domain
220         IF( lk_mpp )   CALL mpp_sum( zsurf    )       ! sum over the global domain
221
222         ! Default GYRE configuration
223         zsumemp = zsumemp / zsurf
224      ELSE
225         ! Benchmark GYRE configuration (to allow the bit to bit comparison between Mpp/Mono case)
226         zsumemp = 0.e0   ;    zsurf = 0.e0
227      ENDIF
228
229      !salinity terms
230      emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)
231      emps(:,:) = emp(:,:)
232
233
234      ! ---------------------------- !
235      !       momentum fluxes        !
236      ! ---------------------------- !
237      ! same wind as in Wico
238      !test date0 : ndate0 = 010203
239      zyear0  =   ndate0 / 10000
240      zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100
241      zday0   =   ndate0 - zyear0 * 10000 - zmonth0 * 100
242      !Calculates nday_year, day since january 1st
243      zday_year0 = (zmonth0-1)*30.+zday0
244
245      !accumulates days of previous months of this year
246      ! day (in hours) since january the 1st
247      ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &  ! incrementation in hour
248         &     - (nyear - 1) * rjjhh * zyydd          !  - nber of hours the precedent years
249      ztimemax = ((5.*30.)+21.)* 24.               ! 21th june     in hours
250      ztimemin = ztimemax + rjjhh * zyydd / 2      ! 21th december in hours
251      !                                            ! NB: rjjhh * zyydd / 4 = 1 seasonal cycle in hours
252
253      ! mean intensity at 0.105 ; srqt(2) because projected with 45deg angle
254      ztau = 0.105 / SQRT( 2. )
255      ! seasonal oscillation intensity
256      ztau_sais = 0.015
257      ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )
258      DO jj = 1, jpj
259         DO ji = 1, jpi
260           ! domain from 15deg to 50deg and 1/2 period along 14deg
261           ! so 5/4 of half period with seasonal cycle
262           utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
263           vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
264         END DO
265      END DO
266
267      ! module of wind stress and wind speed at T-point
268      zcoef = 1. / ( zrhoa * zcdrag ) 
269!CDIR NOVERRCHK
270      DO jj = 2, jpjm1
271!CDIR NOVERRCHK
272         DO ji = fs_2, fs_jpim1   ! vect. opt.
273            ztx = utau(ji-1,jj  ) + utau(ji,jj) 
274            zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
275            zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
276            taum(ji,jj) = zmod
277            wndm(ji,jj) = SQRT( zmod * zcoef )
278         END DO
279      END DO
280      CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
281
282      ! ---------------------------------- !
283      !  control print at first time-step  !
284      ! ---------------------------------- !
285      IF( kt == nit000 .AND. lwp ) THEN
286         WRITE(numout,*)
287         WRITE(numout,*)'sbc_gyre : analytical surface fluxes for GYRE configuration'               
288         WRITE(numout,*)'~~~~~~~~ ' 
289         WRITE(numout,*)'           nyear      = ', nyear
290         WRITE(numout,*)'           nmonth     = ', nmonth
291         WRITE(numout,*)'           nday       = ', nday
292         WRITE(numout,*)'           nday_year  = ', nday_year
293         WRITE(numout,*)'           ztime      = ', ztime
294         WRITE(numout,*)'           ztimemax   = ', ztimemax
295         WRITE(numout,*)'           ztimemin   = ', ztimemin
296         WRITE(numout,*)'           ztimemax1  = ', ztimemax1
297         WRITE(numout,*)'           ztimemin1  = ', ztimemin1
298         WRITE(numout,*)'           ztimemax2  = ', ztimemax2
299         WRITE(numout,*)'           ztimemin2  = ', ztimemin2
300         WRITE(numout,*)'           zyear0     = ', zyear0
301         WRITE(numout,*)'           zmonth0    = ', zmonth0
302         WRITE(numout,*)'           zday0      = ', zday0
303         WRITE(numout,*)'           zday_year0 = ', zday_year0
304         WRITE(numout,*)'           zyydd      = ', zyydd
305         WRITE(numout,*)'           zemp_S     = ', zemp_S
306         WRITE(numout,*)'           zemp_N     = ', zemp_N
307         WRITE(numout,*)'           zemp_sais  = ', zemp_sais
308         WRITE(numout,*)'           zTstar     = ', zTstar
309         WRITE(numout,*)'           zsumemp    = ', zsumemp
310         WRITE(numout,*)'           zsurf      = ', zsurf
311         WRITE(numout,*)'           ztrp       = ', ztrp
312         WRITE(numout,*)'           zconv      = ', zconv
313         WRITE(numout,*)'           ndastp     = ', ndastp
314         WRITE(numout,*)'           adatrj     = ', adatrj
315
316      ENDIF
317      !
318   END SUBROUTINE sbc_gyre
319
320   !!======================================================================
321END MODULE sbcana
Note: See TracBrowser for help on using the repository browser.