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.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.F90 @ 6583

Last change on this file since 6583 was 6583, checked in by flavoni, 8 years ago

add module for sbc gyre in usr_def module , see ticket #1692
File size: 11.1 KB
Line 
1MODULE usrdef
2   !!==============================================================================
3   !!                       ***  MODULE usrdef   ***
4   !! User defined module: used like example to define domain, init, sbc, ...
5   !!==============================================================================
6   !! History :  NEMO ! 2016-03  (S. Flavoni)
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   usr_def_hgr       : initialize the horizontal mesh
11   !!   usr_def_ini       : initial state
12   !!----------------------------------------------------------------------
13   USE step_oce       ! module used in the ocean time stepping module (step.F90)
14   USE mppini         ! shared/distributed memory setting (mpp_init routine)             
15   USE phycst         ! physical constants
16   IMPLICIT NONE
17   PRIVATE
18
19
20   PUBLIC usr_def_hgr
21   PUBLIC usr_def_ini
22
23   !!----------------------------------------------------------------------
24   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
25   !! $Id:$
26   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
27   !!----------------------------------------------------------------------
28CONTAINS
29
30   SUBROUTINE usr_def_hgr( kbench, k_cfg , kff    , pff   , & ! Coriolis parameter  (if domain not on the sphere)
31      &                    pglamt, pglamu, pglamv , pglamf, & ! gridpoints position (required)
32      &                    pgphit, pgphiu, pgphiv , pgphif, & !         
33      &                    pe1t  , pe1u  , pe1v   , pe1f  , & ! scale factors (required)
34      &                    pe2t  , pe2u  , pe2v   , pe2f  , & !
35      &                    ke1e2u_v                       )   ! u- & v-surfaces (if gridsize reduction is used in strait(s))
36      !!----------------------------------------------------------------------
37      !!                  ***  ROUTINE usr_def_hgr  ***
38      !!
39      !! ** Purpose :   user defined mesh and Coriolis parameter
40      !!
41      !! ** Method  :   set all intent(out) argument to a proper value
42      !!
43      !!                Here GYRE configuration :
44      !!          Rectangular mid-latitude domain
45      !!          - with axes rotated by 45 degrees
46      !!          - a constant horizontal resolution of 106 km
47      !!          - on a beta-plane
48      !!
49      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees)
50      !!              - define coriolis parameter at f-point if the domain in not on the sphere (on beta-plane)
51      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters)
52      !!              - define u- & v-surfaces (if gridsize reduction is used in some straits) (in m2)
53      !!----------------------------------------------------------------------
54      INTEGER                 , INTENT(in   ) ::   kbench, k_cfg   ! parameter of namelist for benchmark, and dimension of GYRE
55      INTEGER                 , INTENT(  out) ::   kff             ! =1 Coriolis parameter computed here, =0 otherwise
56      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pff             ! Coriolis factor at f-point                  [1/s]
57      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pglamt, pglamu, pglamv, pglamf !  longitude outputs        [degrees]
58      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgphit, pgphiu, pgphiv, pgphif !  latitude outputs         [degrees]
59      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe1t, pe1u, pe1v, pe1f   !  i-scale factors                      [m]
60      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe2t, pe2u, pe2v, pe2f   !  j-scale factors                      [m]
61      INTEGER                 , INTENT(  out) ::   ke1e2u_v                 ! =1 u- & v-surfaces read here, =0 otherwise
62      !!-------------------------------------------------------------------------------
63      INTEGER  ::   ji, jj               ! dummy loop indices
64      REAL(wp) ::   zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1  , ze1deg, zf0 ! local scalars
65      REAL(wp) ::   zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme      ! local scalars
66      !!-------------------------------------------------------------------------------
67      !
68      IF( nn_timing == 1 )  CALL timing_start('usr_def_hgr')
69      !
70      !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration)
71      !
72      IF(lwp) WRITE(numout,*)
73      IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
74      !
75      ! angle 45deg and ze1=106.e+3 / k_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
76      zlam1 = -85._wp
77      zphi1 =  29._wp
78      ! resolution in meters
79      ze1 = 106000. / REAL( k_cfg , wp )
80      ! benchmark: forced the resolution to be about 100 km
81            IF( kbench /= 0 )   ze1 = 106000._wp   
82      zsin_alpha = - SQRT( 2._wp ) * 0.5_wp
83      zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
84      ze1deg = ze1 / (ra * rad)
85      IF( kbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon
86      !                                                           ! at the right jp_cfg resolution
87      zlam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
88      zphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
89      !   
90      IF( nprint==1 .AND. lwp )   THEN
91         WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
92         WRITE(numout,*) '          ze1deg', ze1deg, 'zlam0', zlam0, 'zphi0', zphi0
93      ENDIF
94      !   
95      DO jj = 1, jpj 
96         DO ji = 1, jpi 
97            zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5 
98            zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5 
99            !   
100            !glamt(i,j) longitude at T-point
101            !gphit(i,j) latitude at T-point 
102            pglamt(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
103            pgphit(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
104            !   
105            !glamu(i,j) longitude at U-point
106            !gphiu(i,j) latitude at U-point
107            pglamu(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
108            pgphiu(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
109            !   
110            !glamv(i,j) longitude at V-point
111            !gphiv(i,j) latitude at V-point
112            pglamv(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
113            pgphiv(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
114            !
115            !glamf(i,j) longitude at F-point
116            !gphif(i,j) latitude at F-point
117            pglamf(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
118            pgphif(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
119         END DO
120      END DO
121      !
122      !       !== Horizontal scale factors ==! (in meters)
123      !                     
124      !                             ! constant grid spacing
125      pe1t(:,:) =  ze1     ;      pe2t(:,:) = ze1
126      pe1u(:,:) =  ze1     ;      pe2u(:,:) = ze1
127      pe1v(:,:) =  ze1     ;      pe2v(:,:) = ze1
128      pe1f(:,:) =  ze1     ;      pe2f(:,:) = ze1
129      !                             ! NO reduction of grid size in some straits
130      ke1e2u_v = 0                  !    ==>> u_ & v_surfaces will be computed in dom_ghr routine
131      !
132      !
133
134      !                      !==  Coriolis parameter  ==!
135      kff = 1                                                           !  indicate not to compute ff afterward
136      !
137      zbeta = 2. * omega * COS( rad * zphi1 ) / ra                      ! beta at latitude zphi1
138      !SF we overwrite zphi0 (south point in latitude) used just above to define pgphif (value of zphi0=15.5190567531966)
139      !SF for computation of coriolis we keep the parameter of Hazeleger, W., and S. S. Drijfhout, JPO 1998.
140      zphi0 = 15._wp                                                    !  latitude of the most southern grid point 
141      zf0   = 2. * omega * SIN( rad * zphi0 )                           !  compute f0 1st point south
142      !
143      pff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south)
144      !
145      IF(lwp) WRITE(numout,*) '                           beta-plane used. beta = ', zbeta, ' 1/(s.m)'
146      !
147      IF( nn_timing == 1 )  CALL timing_stop('usr_def_hgr')
148      !
149   END SUBROUTINE usr_def_hgr
150
151
152   SUBROUTINE usr_def_zgr()           ! Coriolis parameter  (if domain not on the sphere)
153       ! subroutine for vertical grid
154   END SUBROUTINE usr_def_zgr
155
156 
157   SUBROUTINE usr_def_ini( pts )
158      !!----------------------------------------------------------------------
159      !!                   ***  ROUTINE usr_def_ini  ***
160      !!
161      !! ** Purpose :   Initialization of the dynamics and tracers
162      !!                Here GYRE configuration example : (double gyre with rotated domain)
163      !!
164      !! ** Method  : - set temprature field
165      !!              - set salinity   field
166      !!----------------------------------------------------------------------
167      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pts   ! T & S data
168      !
169      INTEGER :: ji, jj, jk  ! dummy loop indices
170      !!----------------------------------------------------------------------
171      !
172      IF(lwp) WRITE(numout,*)
173      IF(lwp) WRITE(numout,*) 'usr_def_ini : analytical definition of initial state '
174      IF(lwp) WRITE(numout,*) '              T and S profiles deduced from LEVITUS '
175      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
176      !
177      DO jk = 1, jpk
178         DO jj = 1, jpj
179            DO ji = 1, jpi
180               pts(ji,jj,jk,jp_tem) =  (  (  16. - 12. * TANH( (gdept_n(ji,jj,jk) - 400) / 700 ) )   &
181                    &           * (-TANH( (500-gdept_n(ji,jj,jk)) / 150 ) + 1) / 2               &
182                    &           + ( 15. * ( 1. - TANH( (gdept_n(ji,jj,jk)-50.) / 1500.) )        &
183                    &           - 1.4 * TANH((gdept_n(ji,jj,jk)-100.) / 100.)                    &
184                    &           + 7.  * (1500. - gdept_n(ji,jj,jk)) / 1500.)                     &
185                    &           * (-TANH( (gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2  ) * tmask(ji,jj,jk)
186
187               pts(ji,jj,jk,jp_sal) =  (  (  36.25 - 1.13 * TANH( (gdept_n(ji,jj,jk) - 305) / 460 ) )  &
188                    &         * (-TANH((500 - gdept_n(ji,jj,jk)) / 150) + 1) / 2              &
189                    &         + ( 35.55 + 1.25 * (5000. - gdept_n(ji,jj,jk)) / 5000.          &
190                    &         - 1.62 * TANH( (gdept_n(ji,jj,jk) - 60.  ) / 650. )             &
191                    &         + 0.2  * TANH( (gdept_n(ji,jj,jk) - 35.  ) / 100. )             &
192                    &         + 0.2  * TANH( (gdept_n(ji,jj,jk) - 1000.) / 5000.) )           &
193                    &         * (-TANH((gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2  ) * tmask(ji,jj,jk)
194            END DO
195         END DO
196      END DO
197      !   
198   END SUBROUTINE usr_def_ini
199
200   !!======================================================================
201END MODULE usrdef
Note: See TracBrowser for help on using the repository browser.