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/NEMO/OPA_SRC/SBC – NEMO

source: trunk/NEMO/OPA_SRC/SBC/sbcana.F90 @ 702

Last change on this file since 702 was 702, checked in by smasson, 17 years ago

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