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 NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/USR – NEMO

source: NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/USR/usrdef_sbc.F90 @ 15432

Last change on this file since 15432 was 15432, checked in by jcastill, 12 months ago

Reverting the changes in USR as they are now included in the CO9 branch

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