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.
usrdef_sbc.F90 in branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC – NEMO

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef_sbc.F90 @ 6595

Last change on this file since 6595 was 6595, checked in by gm, 8 years ago

#1692 - branch SIMPLIF_2_usrdef: add namusr_def, domain size read in file, interface with arguments for hgr_read and usr_def_hgr

File size: 11.5 KB
Line 
1MODULE usrdef_sbc
2   !!======================================================================
3   !!                       ***  MODULE  usrdef_sbc  ***
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   !!   usr_def_sbc    : user defined surface bounday conditions
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   !
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   usr_def_sbc    ! routine called in sbcmod module
27
28   !! * Substitutions
29#  include "vectopt_loop_substitute.h90"
30   !!----------------------------------------------------------------------
31   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
32   !! $Id:$
33   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
34   !!----------------------------------------------------------------------
35CONTAINS
36
37   SUBROUTINE usr_def_sbc( kt )
38      !!---------------------------------------------------------------------
39      !!                    ***  ROUTINE usr_def_sbc  ***
40      !!             
41      !! ** Purpose :   provide at each time-step the GYRE surface boundary
42      !!              condition, i.e. the momentum, heat and freshwater fluxes.
43      !!
44      !! ** Method  :   analytical seasonal cycle for GYRE configuration.
45      !!                CAUTION : never mask the surface stress field !
46      !!
47      !! ** Action  : - set the ocean surface boundary condition, i.e.   
48      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
49      !!
50      !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000.
51      !!----------------------------------------------------------------------
52      INTEGER, INTENT(in) ::   kt   ! ocean time step
53      !!
54      INTEGER  ::   ji, jj                 ! dummy loop indices
55      INTEGER  ::   zyear0                 ! initial year
56      INTEGER  ::   zmonth0                ! initial month
57      INTEGER  ::   zday0                  ! initial day
58      INTEGER  ::   zday_year0             ! initial day since january 1st
59      REAL(wp) ::   ztau     , ztau_sais   ! wind intensity and of the seasonal cycle
60      REAL(wp) ::   ztime                  ! time in hour
61      REAL(wp) ::   ztimemax , ztimemin    ! 21th June, and 21th decem. if date0 = 1st january
62      REAL(wp) ::   ztimemax1, ztimemin1   ! 21th June, and 21th decem. if date0 = 1st january
63      REAL(wp) ::   ztimemax2, ztimemin2   ! 21th June, and 21th decem. if date0 = 1st january
64      REAL(wp) ::   ztaun                  ! intensity
65      REAL(wp) ::   zemp_s, zemp_n, zemp_sais, ztstar
66      REAL(wp) ::   zcos_sais1, zcos_sais2, ztrp, zconv, t_star
67      REAL(wp) ::   zsumemp, zsurf
68      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
69      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
70      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
71      REAL(wp) ::   zyydd                 ! number of days in one year
72      !!---------------------------------------------------------------------
73      zyydd = REAL(nyear_len(1),wp)
74
75      ! ---------------------------- !
76      !  heat and freshwater fluxes  !
77      ! ---------------------------- !
78      !same temperature, E-P as in HAZELEGER 2000
79
80      zyear0     =   ndate0 / 10000                             ! initial year
81      zmonth0    = ( ndate0 - zyear0 * 10000 ) / 100            ! initial month
82      zday0      =   ndate0 - zyear0 * 10000 - zmonth0 * 100    ! initial day betwen 1 and 30
83      zday_year0 = ( zmonth0 - 1 ) * 30.+zday0                  ! initial day betwen 1 and 360
84
85      ! current day (in hours) since january the 1st of the current year
86      ztime = REAL( kt ) * rdt / (rmmss * rhhmm)   &       !  total incrementation (in hours)
87         &      - (nyear  - 1) * rjjhh * zyydd             !  minus years since beginning of experiment (in hours)
88
89      ztimemax1 = ((5.*30.)+21.)* 24.                      ! 21th june     at 24h in hours
90      ztimemin1 = ztimemax1 + rjjhh * zyydd / 2            ! 21th december        in hours
91      ztimemax2 = ((6.*30.)+21.)* 24.                      ! 21th july     at 24h in hours
92      ztimemin2 = ztimemax2 - rjjhh * zyydd / 2            ! 21th january         in hours
93      !                                                    ! NB: rjjhh * zyydd / 4 = one seasonal cycle in hours
94
95      ! amplitudes
96      zemp_S    = 0.7       ! intensity of COS in the South
97      zemp_N    = 0.8       ! intensity of COS in the North
98      zemp_sais = 0.1
99      zTstar    = 28.3      ! intemsity from 28.3 a -5 deg
100
101      ! 1/2 period between 21th June and 21th December and between 21th July and 21th January
102      zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi ) 
103      zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi )
104
105      ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K)
106      zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s
107      DO jj = 1, jpj
108         DO ji = 1, jpi
109            ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
110            ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
111            ! 64.5 in summer, 42.5 in winter
112            t_star = zTstar * ( 1. + 1. / 50. * zcos_sais2 )                &
113               &                    * COS( rpi * (gphit(ji,jj) - 5.)               &
114               &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) )
115            ! 23.5 deg : tropics
116            qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
117            qns (ji,jj) = ztrp * ( tsb(ji,jj,1,jp_tem) - t_star ) - qsr(ji,jj)
118            IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg
119               emp  (ji,jj) =   zemp_S * zconv   &
120                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  &
121                  &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
122            ELSE
123               emp (ji,jj) =  - zemp_N * zconv   &
124                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  &
125                  &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
126            ENDIF
127         END DO
128      END DO
129
130!!gm  ===>>> The following lines are to be de-commented
131!!    Commented here as they change the GYRE results
132!#if !defined key_mpp_rep
133!      CALL ctl_stop( 'usr_def_sbc: key_mpp_rep is required to run GYRE configuration')
134!#endif
135!!gm
136
137      zsumemp = GLOB_SUM( emp  (:,:)   ) 
138      zsurf   = GLOB_SUM( tmask(:,:,1) ) 
139      zsumemp = zsumemp / zsurf         ! Default GYRE configuration
140
141      ! freshwater (mass flux) and update of qns with heat content of emp
142      emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average)
143      sfx (:,:) = 0.0_wp                                   ! no salt flux
144      qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST
145
146
147      ! ---------------------------- !
148      !       momentum fluxes        !
149      ! ---------------------------- !
150      ! same wind as in Wico
151      !test date0 : ndate0 = 010203
152      zyear0  =   ndate0 / 10000
153      zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100
154      zday0   =   ndate0 - zyear0 * 10000 - zmonth0 * 100
155      !Calculates nday_year, day since january 1st
156      zday_year0 = (zmonth0-1)*30.+zday0
157
158      !accumulates days of previous months of this year
159      ! day (in hours) since january the 1st
160      ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &  ! incrementation in hour
161         &     - (nyear - 1) * rjjhh * zyydd          !  - nber of hours the precedent years
162      ztimemax = ((5.*30.)+21.)* 24.               ! 21th june     in hours
163      ztimemin = ztimemax + rjjhh * zyydd / 2      ! 21th december in hours
164      !                                            ! NB: rjjhh * zyydd / 4 = 1 seasonal cycle in hours
165
166      ! mean intensity at 0.105 ; srqt(2) because projected with 45deg angle
167      ztau = 0.105 / SQRT( 2. )
168      ! seasonal oscillation intensity
169      ztau_sais = 0.015
170      ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )
171      DO jj = 1, jpj
172         DO ji = 1, jpi
173           ! domain from 15deg to 50deg and 1/2 period along 14deg
174           ! so 5/4 of half period with seasonal cycle
175           utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
176           vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
177         END DO
178      END DO
179
180      ! module of wind stress and wind speed at T-point
181      zcoef = 1. / ( zrhoa * zcdrag ) 
182      DO jj = 2, jpjm1
183         DO ji = fs_2, fs_jpim1   ! vect. opt.
184            ztx = utau(ji-1,jj  ) + utau(ji,jj) 
185            zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
186            zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
187            taum(ji,jj) = zmod
188            wndm(ji,jj) = SQRT( zmod * zcoef )
189         END DO
190      END DO
191      CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
192
193      ! ---------------------------------- !
194      !  control print at first time-step  !
195      ! ---------------------------------- !
196      IF( kt == nit000 .AND. lwp ) THEN
197         WRITE(numout,*)
198         WRITE(numout,*)'usr_def_sbc : analytical surface fluxes for GYRE configuration'               
199         WRITE(numout,*)'~~~~~~~~~~~ ' 
200         WRITE(numout,*)'           nyear      = ', nyear
201         WRITE(numout,*)'           nmonth     = ', nmonth
202         WRITE(numout,*)'           nday       = ', nday
203         WRITE(numout,*)'           nday_year  = ', nday_year
204         WRITE(numout,*)'           ztime      = ', ztime
205         WRITE(numout,*)'           ztimemax   = ', ztimemax
206         WRITE(numout,*)'           ztimemin   = ', ztimemin
207         WRITE(numout,*)'           ztimemax1  = ', ztimemax1
208         WRITE(numout,*)'           ztimemin1  = ', ztimemin1
209         WRITE(numout,*)'           ztimemax2  = ', ztimemax2
210         WRITE(numout,*)'           ztimemin2  = ', ztimemin2
211         WRITE(numout,*)'           zyear0     = ', zyear0
212         WRITE(numout,*)'           zmonth0    = ', zmonth0
213         WRITE(numout,*)'           zday0      = ', zday0
214         WRITE(numout,*)'           zday_year0 = ', zday_year0
215         WRITE(numout,*)'           zyydd      = ', zyydd
216         WRITE(numout,*)'           zemp_S     = ', zemp_S
217         WRITE(numout,*)'           zemp_N     = ', zemp_N
218         WRITE(numout,*)'           zemp_sais  = ', zemp_sais
219         WRITE(numout,*)'           zTstar     = ', zTstar
220         WRITE(numout,*)'           zsumemp    = ', zsumemp
221         WRITE(numout,*)'           zsurf      = ', zsurf
222         WRITE(numout,*)'           ztrp       = ', ztrp
223         WRITE(numout,*)'           zconv      = ', zconv
224         WRITE(numout,*)'           ndastp     = ', ndastp
225         WRITE(numout,*)'           adatrj     = ', adatrj
226      ENDIF
227      !
228   END SUBROUTINE usr_def_sbc
229
230   !!======================================================================
231END MODULE usrdef_sbc
Note: See TracBrowser for help on using the repository browser.