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/2019/dev_r11943_MERGE_2019/src/OCE/USR – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/USR/usrdef_sbc.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 11.8 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#  include "do_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
37   !! $Id$
38   !! Software governed by the CeCILL license (see ./LICENSE)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE usrdef_sbc_oce( kt, Kbb )
43      !!---------------------------------------------------------------------
44      !!                    ***  ROUTINE usrdef_sbc  ***
45      !!             
46      !! ** Purpose :   provide at each time-step the GYRE surface boundary
47      !!              condition, i.e. the momentum, heat and freshwater fluxes.
48      !!
49      !! ** Method  :   analytical seasonal cycle for GYRE configuration.
50      !!                CAUTION : never mask the surface stress field !
51      !!
52      !! ** Action  : - set the ocean surface boundary condition, i.e.   
53      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
54      !!
55      !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000.
56      !!----------------------------------------------------------------------
57      INTEGER, INTENT(in) ::   kt   ! ocean time step
58      INTEGER, INTENT(in) ::   Kbb  ! ocean time index
59      !!
60      INTEGER  ::   ji, jj                 ! dummy loop indices
61      INTEGER  ::   zyear0                 ! initial year
62      INTEGER  ::   zmonth0                ! initial month
63      INTEGER  ::   zday0                  ! initial day
64      INTEGER  ::   zday_year0             ! initial day since january 1st
65      REAL(wp) ::   ztau     , ztau_sais   ! wind intensity and of the seasonal cycle
66      REAL(wp) ::   ztime                  ! time in hour
67      REAL(wp) ::   ztimemax , ztimemin    ! 21th June, and 21th decem. if date0 = 1st january
68      REAL(wp) ::   ztimemax1, ztimemin1   ! 21th June, and 21th decem. if date0 = 1st january
69      REAL(wp) ::   ztimemax2, ztimemin2   ! 21th June, and 21th decem. if date0 = 1st january
70      REAL(wp) ::   ztaun                  ! intensity
71      REAL(wp) ::   zemp_s, zemp_n, zemp_sais, ztstar
72      REAL(wp) ::   zcos_sais1, zcos_sais2, ztrp, zconv, t_star
73      REAL(wp) ::   zsumemp, zsurf
74      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
75      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
76      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
77      REAL(wp) ::   zyydd                 ! number of days in one year
78      !!---------------------------------------------------------------------
79      zyydd = REAL(nyear_len(1),wp)
80
81      ! ---------------------------- !
82      !  heat and freshwater fluxes  !
83      ! ---------------------------- !
84      !same temperature, E-P as in HAZELEGER 2000
85
86      zyear0     =   ndate0 / 10000                             ! initial year
87      zmonth0    = ( ndate0 - zyear0 * 10000 ) / 100            ! initial month
88      zday0      =   ndate0 - zyear0 * 10000 - zmonth0 * 100    ! initial day betwen 1 and 30
89      zday_year0 = ( zmonth0 - 1 ) * 30.+zday0                  ! initial day betwen 1 and 360
90
91      ! current day (in hours) since january the 1st of the current year
92      ztime = REAL( kt ) * rdt / (rmmss * rhhmm)   &       !  total incrementation (in hours)
93         &      - (nyear  - 1) * rjjhh * zyydd             !  minus years since beginning of experiment (in hours)
94
95      ztimemax1 = ((5.*30.)+21.)* 24.                      ! 21th june     at 24h in hours
96      ztimemin1 = ztimemax1 + rjjhh * zyydd / 2            ! 21th december        in hours
97      ztimemax2 = ((6.*30.)+21.)* 24.                      ! 21th july     at 24h in hours
98      ztimemin2 = ztimemax2 - rjjhh * zyydd / 2            ! 21th january         in hours
99      !                                                    ! NB: rjjhh * zyydd / 4 = one seasonal cycle in hours
100
101      ! amplitudes
102      zemp_S    = 0.7       ! intensity of COS in the South
103      zemp_N    = 0.8       ! intensity of COS in the North
104      zemp_sais = 0.1
105      zTstar    = 28.3      ! intemsity from 28.3 a -5 deg
106
107      ! 1/2 period between 21th June and 21th December and between 21th July and 21th January
108      zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi ) 
109      zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi )
110
111      ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K)
112      zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s
113      DO_2D_11_11
114         ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
115         ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
116         ! 64.5 in summer, 42.5 in winter
117         t_star = zTstar * ( 1. + 1. / 50. * zcos_sais2 )                &
118            &                    * COS( rpi * (gphit(ji,jj) - 5.)               &
119            &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) )
120         ! 23.5 deg : tropics
121         qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
122         qns (ji,jj) = ztrp * ( ts(ji,jj,1,jp_tem,Kbb) - t_star ) - qsr(ji,jj)
123         IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg
124            emp  (ji,jj) =   zemp_S * zconv   &
125               &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  &
126               &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
127         ELSE
128            emp (ji,jj) =  - zemp_N * zconv   &
129               &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  &
130               &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
131         ENDIF
132      END_2D
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_2D_11_11
169        ! domain from 15deg to 50deg and 1/2 period along 14deg
170        ! so 5/4 of half period with seasonal cycle
171        utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
172        vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
173      END_2D
174
175      ! module of wind stress and wind speed at T-point
176      zcoef = 1. / ( zrhoa * zcdrag ) 
177      DO_2D_00_00
178         ztx = utau(ji-1,jj  ) + utau(ji,jj) 
179         zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
180         zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
181         taum(ji,jj) = zmod
182         wndm(ji,jj) = SQRT( zmod * zcoef )
183      END_2D
184      CALL lbc_lnk_multi( 'usrdef_sbc', taum(:,:), 'T', 1. , wndm(:,:), 'T', 1. )
185
186      ! ---------------------------------- !
187      !  control print at first time-step  !
188      ! ---------------------------------- !
189      IF( kt == nit000 .AND. lwp ) THEN
190         WRITE(numout,*)
191         WRITE(numout,*)'usrdef_sbc_oce : analytical surface fluxes for GYRE configuration'               
192         WRITE(numout,*)'~~~~~~~~~~~ ' 
193         WRITE(numout,*)'           nyear      = ', nyear
194         WRITE(numout,*)'           nmonth     = ', nmonth
195         WRITE(numout,*)'           nday       = ', nday
196         WRITE(numout,*)'           nday_year  = ', nday_year
197         WRITE(numout,*)'           ztime      = ', ztime
198         WRITE(numout,*)'           ztimemax   = ', ztimemax
199         WRITE(numout,*)'           ztimemin   = ', ztimemin
200         WRITE(numout,*)'           ztimemax1  = ', ztimemax1
201         WRITE(numout,*)'           ztimemin1  = ', ztimemin1
202         WRITE(numout,*)'           ztimemax2  = ', ztimemax2
203         WRITE(numout,*)'           ztimemin2  = ', ztimemin2
204         WRITE(numout,*)'           zyear0     = ', zyear0
205         WRITE(numout,*)'           zmonth0    = ', zmonth0
206         WRITE(numout,*)'           zday0      = ', zday0
207         WRITE(numout,*)'           zday_year0 = ', zday_year0
208         WRITE(numout,*)'           zyydd      = ', zyydd
209         WRITE(numout,*)'           zemp_S     = ', zemp_S
210         WRITE(numout,*)'           zemp_N     = ', zemp_N
211         WRITE(numout,*)'           zemp_sais  = ', zemp_sais
212         WRITE(numout,*)'           zTstar     = ', zTstar
213         WRITE(numout,*)'           zsumemp    = ', zsumemp
214         WRITE(numout,*)'           zsurf      = ', zsurf
215         WRITE(numout,*)'           ztrp       = ', ztrp
216         WRITE(numout,*)'           zconv      = ', zconv
217         WRITE(numout,*)'           ndastp     = ', ndastp
218         WRITE(numout,*)'           adatrj     = ', adatrj
219      ENDIF
220      !
221   END SUBROUTINE usrdef_sbc_oce
222
223
224   SUBROUTINE usrdef_sbc_ice_tau( kt )
225      INTEGER, INTENT(in) ::   kt   ! ocean time step
226   END SUBROUTINE usrdef_sbc_ice_tau
227
228
229   SUBROUTINE usrdef_sbc_ice_flx( kt, phs, phi )
230      INTEGER, INTENT(in) ::   kt   ! ocean time step
231      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
232      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
233   END SUBROUTINE usrdef_sbc_ice_flx
234
235   !!======================================================================
236END MODULE usrdef_sbc
Note: See TracBrowser for help on using the repository browser.