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 branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90 @ 5550

Last change on this file since 5550 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.

  • Property svn:keywords set to Id
File size: 16.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   USE lib_fortran
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  ! nb of time-step during which the surface stress
31   !                        ! increase from 0 to its nominal value
32   REAL(wp) ::   rn_utau0   ! constant wind stress value in i-direction
33   REAL(wp) ::   rn_vtau0   ! constant wind stress value in j-direction
34   REAL(wp) ::   rn_qns0    ! non solar heat flux
35   REAL(wp) ::   rn_qsr0    !     solar heat flux
36   REAL(wp) ::   rn_emp0    ! net freshwater flux
37   
38   !! * Substitutions
39#  include "domzgr_substitute.h90"
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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      !!
60      !! ** Action  : - set the ocean surface boundary condition, i.e. 
61      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
62      !!----------------------------------------------------------------------
63      INTEGER, INTENT(in) ::   kt   ! ocean time step
64      !
65      INTEGER  ::   ios                   ! local integer
66      REAL(wp) ::   zrhoa  = 1.22_wp      ! air density kg/m3
67      REAL(wp) ::   zcdrag = 1.5e-3_wp    ! drag coefficient
68      REAL(wp) ::   zfact, ztx            ! local scalars
69      REAL(wp) ::   zcoef, zty, zmod      !   -      -
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_ref )              ! Namelist namsbc_ana in reference namelist : Analytical surface fluxes
77         READ  ( numnam_ref, namsbc_ana, IOSTAT = ios, ERR = 901)
78901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_ana in reference namelist', lwp )
79
80         REWIND( numnam_cfg )              ! Namelist namsbc_ana in configuration namelist : Analytical surface fluxes
81         READ  ( numnam_cfg, namsbc_ana, IOSTAT = ios, ERR = 902 )
82902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_ana in configuration namelist', lwp )
83         IF(lwm) WRITE ( numond, namsbc_ana )
84         !
85         IF(lwp) WRITE(numout,*)' '
86         IF(lwp) WRITE(numout,*)' sbc_ana : Constant surface fluxes read in namsbc_ana namelist'
87         IF(lwp) WRITE(numout,*)' ~~~~~~~ '
88         IF(lwp) WRITE(numout,*)'              spin up of the stress  nn_tau000 = ', nn_tau000, ' time-steps'
89         IF(lwp) WRITE(numout,*)'              constant i-stress      rn_utau0  = ', rn_utau0 , ' N/m2'
90         IF(lwp) WRITE(numout,*)'              constant j-stress      rn_vtau0  = ', rn_vtau0 , ' N/m2'
91         IF(lwp) WRITE(numout,*)'              non solar heat flux    rn_qns0   = ', rn_qns0  , ' W/m2'
92         IF(lwp) WRITE(numout,*)'              solar heat flux        rn_qsr0   = ', rn_qsr0  , ' W/m2'
93         IF(lwp) WRITE(numout,*)'              net heat flux          rn_emp0   = ', rn_emp0  , ' Kg/m2/s'
94         !
95         nn_tau000 = MAX( nn_tau000, 1 )     ! must be >= 1
96         !
97         utau(:,:) = rn_utau0
98         vtau(:,:) = rn_vtau0
99         taum(:,:) = SQRT ( rn_utau0 * rn_utau0 + rn_vtau0 * rn_vtau0 )
100         wndm(:,:) = SQRT ( taum(1,1) /  ( zrhoa * zcdrag ) )
101         !
102         emp (:,:) = rn_emp0
103         sfx (:,:) = 0.0_wp
104         qns (:,:) = rn_qns0 - emp(:,:) * sst_m(:,:) * rcp      ! including heat content associated with mass flux at SST
105         qsr (:,:) = rn_qsr0
106         !         
107      ENDIF
108
109      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
110         !
111         IF( kt <= nn_tau000 ) THEN       ! Increase the stress to its nominal value
112            !                             ! during the first nn_tau000 time-steps
113            zfact = 0.5 * (  1. - COS( rpi * REAL( kt, wp ) / REAL( nn_tau000, wp ) )  )
114            zcoef = 1. / ( zrhoa * zcdrag ) 
115            ztx   = zfact * rn_utau0
116            zty   = zfact * rn_vtau0
117            zmod  = SQRT( ztx * ztx + zty * zty )
118            utau(:,:) = ztx
119            vtau(:,:) = zty
120            taum(:,:) = zmod
121            zmod = SQRT( zmod * zcoef )
122            wndm(:,:) = zmod
123         ENDIF
124         !                                ! update heat and fresh water fluxes
125         !                                ! as they may have been changed by sbcssr module
126         emp (:,:) = rn_emp0              ! NB: qns changes with SST if emp /= 0
127         sfx (:,:) = 0._wp
128         qns (:,:) = rn_qns0 - emp(:,:) * sst_m(:,:) * rcp
129         qsr (:,:) = rn_qsr0
130         !
131      ENDIF
132      !
133   END SUBROUTINE sbc_ana
134
135
136   SUBROUTINE sbc_gyre( kt )
137      !!---------------------------------------------------------------------
138      !!                    ***  ROUTINE sbc_ana ***
139      !!             
140      !! ** Purpose :   provide at each time-step the GYRE surface boundary
141      !!              condition, i.e. the momentum, heat and freshwater fluxes.
142      !!
143      !! ** Method  :   analytical seasonal cycle for GYRE configuration.
144      !!                CAUTION : never mask the surface stress field !
145      !!
146      !! ** Action  : - set the ocean surface boundary condition, i.e.   
147      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
148      !!
149      !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000.
150      !!----------------------------------------------------------------------
151      INTEGER, INTENT(in) ::   kt          ! ocean time step
152      !!
153      INTEGER  ::   ji, jj                 ! dummy loop indices
154      INTEGER  ::   zyear0                 ! initial year
155      INTEGER  ::   zmonth0                ! initial month
156      INTEGER  ::   zday0                  ! initial day
157      INTEGER  ::   zday_year0             ! initial day since january 1st
158      REAL(wp) ::   ztau     , ztau_sais   ! wind intensity and of the seasonal cycle
159      REAL(wp) ::   ztime                  ! time in hour
160      REAL(wp) ::   ztimemax , ztimemin    ! 21th June, and 21th decem. if date0 = 1st january
161      REAL(wp) ::   ztimemax1, ztimemin1   ! 21th June, and 21th decem. if date0 = 1st january
162      REAL(wp) ::   ztimemax2, ztimemin2   ! 21th June, and 21th decem. if date0 = 1st january
163      REAL(wp) ::   ztaun                  ! intensity
164      REAL(wp) ::   zemp_s, zemp_n, zemp_sais, ztstar
165      REAL(wp) ::   zcos_sais1, zcos_sais2, ztrp, zconv, t_star
166      REAL(wp) ::   zsumemp, zsurf
167      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
168      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
169      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
170      REAL(wp) ::   zyydd                 ! number of days in one year
171      !!---------------------------------------------------------------------
172      zyydd = REAL(nyear_len(1),wp)
173
174      ! ---------------------------- !
175      !  heat and freshwater fluxes  !
176      ! ---------------------------- !
177      !same temperature, E-P as in HAZELEGER 2000
178
179      zyear0     =   ndate0 / 10000                             ! initial year
180      zmonth0    = ( ndate0 - zyear0 * 10000 ) / 100            ! initial month
181      zday0      =   ndate0 - zyear0 * 10000 - zmonth0 * 100    ! initial day betwen 1 and 30
182      zday_year0 = ( zmonth0 - 1 ) * 30.+zday0                  ! initial day betwen 1 and 360
183
184      ! current day (in hours) since january the 1st of the current year
185      ztime = REAL( kt ) * rdt / (rmmss * rhhmm)   &       !  total incrementation (in hours)
186         &      - (nyear  - 1) * rjjhh * zyydd             !  minus years since beginning of experiment (in hours)
187
188      ztimemax1 = ((5.*30.)+21.)* 24.                      ! 21th june     at 24h in hours
189      ztimemin1 = ztimemax1 + rjjhh * zyydd / 2            ! 21th december        in hours
190      ztimemax2 = ((6.*30.)+21.)* 24.                      ! 21th july     at 24h in hours
191      ztimemin2 = ztimemax2 - rjjhh * zyydd / 2            ! 21th january         in hours
192      !                                                    ! NB: rjjhh * zyydd / 4 = one seasonal cycle in hours
193
194      ! amplitudes
195      zemp_S    = 0.7       ! intensity of COS in the South
196      zemp_N    = 0.8       ! intensity of COS in the North
197      zemp_sais = 0.1
198      zTstar    = 28.3      ! intemsity from 28.3 a -5 deg
199
200      ! 1/2 period between 21th June and 21th December and between 21th July and 21th January
201      zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi ) 
202      zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi )
203
204      ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K)
205      zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s
206      DO jj = 1, jpj
207         DO ji = 1, jpi
208            ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
209            ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
210            ! 64.5 in summer, 42.5 in winter
211            t_star = zTstar * ( 1 + 1. / 50. * zcos_sais2 )                &
212               &                    * COS( rpi * (gphit(ji,jj) - 5.)               &
213               &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) )
214            ! 23.5 deg : tropics
215            qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
216            qns (ji,jj) = ztrp * ( tsb(ji,jj,1,jp_tem) - t_star ) - qsr(ji,jj)
217            IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg
218               emp  (ji,jj) =   zemp_S * zconv   &
219                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  &
220                  &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
221            ELSE
222               emp (ji,jj) =  - zemp_N * zconv   &
223                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  &
224                  &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
225            ENDIF
226         END DO
227      END DO
228
229      ! Compute the emp flux such as its integration on the whole domain at each time is zero
230      IF( nbench /= 1 ) THEN
231         zsumemp = GLOB_SUM( emp(:,:) ) 
232         zsurf   = GLOB_SUM( tmask(:,:,1) ) 
233         ! Default GYRE configuration
234         zsumemp = zsumemp / zsurf
235      ELSE
236         ! Benchmark GYRE configuration (to allow the bit to bit comparison between Mpp/Mono case)
237         zsumemp = 0.e0   ;    zsurf = 0.e0
238      ENDIF
239
240      ! freshwater (mass flux) and update of qns with heat content of emp
241      emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average)
242      sfx (:,:) = 0.0_wp                                   ! no salt flux
243      qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST
244
245
246      ! ---------------------------- !
247      !       momentum fluxes        !
248      ! ---------------------------- !
249      ! same wind as in Wico
250      !test date0 : ndate0 = 010203
251      zyear0  =   ndate0 / 10000
252      zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100
253      zday0   =   ndate0 - zyear0 * 10000 - zmonth0 * 100
254      !Calculates nday_year, day since january 1st
255      zday_year0 = (zmonth0-1)*30.+zday0
256
257      !accumulates days of previous months of this year
258      ! day (in hours) since january the 1st
259      ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &  ! incrementation in hour
260         &     - (nyear - 1) * rjjhh * zyydd          !  - nber of hours the precedent years
261      ztimemax = ((5.*30.)+21.)* 24.               ! 21th june     in hours
262      ztimemin = ztimemax + rjjhh * zyydd / 2      ! 21th december in hours
263      !                                            ! NB: rjjhh * zyydd / 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      ! module of wind stress and wind speed at T-point
280      zcoef = 1. / ( zrhoa * zcdrag ) 
281!CDIR NOVERRCHK
282      DO jj = 2, jpjm1
283!CDIR NOVERRCHK
284         DO ji = fs_2, fs_jpim1   ! vect. opt.
285            ztx = utau(ji-1,jj  ) + utau(ji,jj) 
286            zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
287            zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
288            taum(ji,jj) = zmod
289            wndm(ji,jj) = SQRT( zmod * zcoef )
290         END DO
291      END DO
292      CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
293
294      ! ---------------------------------- !
295      !  control print at first time-step  !
296      ! ---------------------------------- !
297      IF( kt == nit000 .AND. lwp ) THEN
298         WRITE(numout,*)
299         WRITE(numout,*)'sbc_gyre : analytical surface fluxes for GYRE configuration'               
300         WRITE(numout,*)'~~~~~~~~ ' 
301         WRITE(numout,*)'           nyear      = ', nyear
302         WRITE(numout,*)'           nmonth     = ', nmonth
303         WRITE(numout,*)'           nday       = ', nday
304         WRITE(numout,*)'           nday_year  = ', nday_year
305         WRITE(numout,*)'           ztime      = ', ztime
306         WRITE(numout,*)'           ztimemax   = ', ztimemax
307         WRITE(numout,*)'           ztimemin   = ', ztimemin
308         WRITE(numout,*)'           ztimemax1  = ', ztimemax1
309         WRITE(numout,*)'           ztimemin1  = ', ztimemin1
310         WRITE(numout,*)'           ztimemax2  = ', ztimemax2
311         WRITE(numout,*)'           ztimemin2  = ', ztimemin2
312         WRITE(numout,*)'           zyear0     = ', zyear0
313         WRITE(numout,*)'           zmonth0    = ', zmonth0
314         WRITE(numout,*)'           zday0      = ', zday0
315         WRITE(numout,*)'           zday_year0 = ', zday_year0
316         WRITE(numout,*)'           zyydd      = ', zyydd
317         WRITE(numout,*)'           zemp_S     = ', zemp_S
318         WRITE(numout,*)'           zemp_N     = ', zemp_N
319         WRITE(numout,*)'           zemp_sais  = ', zemp_sais
320         WRITE(numout,*)'           zTstar     = ', zTstar
321         WRITE(numout,*)'           zsumemp    = ', zsumemp
322         WRITE(numout,*)'           zsurf      = ', zsurf
323         WRITE(numout,*)'           ztrp       = ', ztrp
324         WRITE(numout,*)'           zconv      = ', zconv
325         WRITE(numout,*)'           ndastp     = ', ndastp
326         WRITE(numout,*)'           adatrj     = ', adatrj
327      ENDIF
328      !
329   END SUBROUTINE sbc_gyre
330
331   !!======================================================================
332END MODULE sbcana
Note: See TracBrowser for help on using the repository browser.