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 @ 6596

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

#1692 - branch SIMPLIF_2_usrdef: remove from namcfg and namdom many obsolete variables ; remove izoom/jzoom option
File size: 14.4 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_nam       : read user defined namelist and set global domain size
11   !!   usr_def_hgr       : initialize the horizontal mesh
12   !!   usr_def_ini       : initial state
13   !!----------------------------------------------------------------------
14   USE dom_oce  , ONLY: nimpp, njmpp       ! ocean space and time domain
15   USE par_oce        ! ocean space and time domain
16   USE phycst         ! physical constants
17   !
18   USE in_out_manager ! I/O manager
19   USE lib_mpp        ! MPP library
20   USE timing         ! Timing
21   
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   usr_def_nam   ! called in nemogcm.F90 module
26   PUBLIC   usr_def_hgr   ! called in domhgr.F90  module
27   PUBLIC   usr_def_ini   ! called in istate.F90  module
28
29   !                      !!* namusr_def namelist *!!
30   LOGICAL ::   ln_bench   ! =T benchmark test with gyre: the gridsize is constant (no need to adjust timestep or viscosity)
31   INTEGER ::   nn_GYRE    ! 1/nn_GYRE = the resolution chosen in degrees and thus defining the horizontal domain size
32
33   !!----------------------------------------------------------------------
34   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
35   !! $Id:$
36   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE usr_def_nam( ldtxt, ldnam, kpi, kpj, kpk )
41      !!----------------------------------------------------------------------
42      !!                     ***  ROUTINE dom_nam  ***
43      !!                   
44      !! ** Purpose :   read user defined namelist and define the domain size
45      !!
46      !! ** Method  :   read in namusr_def containing all the user specific namelist parameter
47      !!
48      !!                Here GYRE configuration
49      !!
50      !! ** input   : - namusr_def namelist found in namelist_cfg
51      !!----------------------------------------------------------------------
52      CHARACTER(len=*), DIMENSION(:), INTENT(out) ::   ldtxt, ldnam    ! stored print information
53      INTEGER                       , INTENT(out) ::   kpi, kpj, kpk   ! global domain sizes
54      !
55      INTEGER ::   ios, ii   ! Local integer
56      !!
57      NAMELIST/namusr_def/ nn_GYRE, ln_bench, jpkglo
58      !!----------------------------------------------------------------------
59      !
60      ii = 1
61      !
62      REWIND( numnam_cfg )          ! Namelist namusr_def (exist in namelist_cfg only)
63      READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 )
64902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist', .TRUE. )
65      !
66!!gm  This does not work...  I don't know how to write namusr_def in "output.namelist.dyn"
67      WRITE( ldnam(ii), namusr_def )
68      !
69      kpi = 30 * nn_GYRE + 2        ! Global Domain size
70      kpj = 20 * nn_GYRE + 2
71      kpk = jpkglo
72      !
73      !                             ! control print
74      WRITE(ldtxt(ii),*)                                                                              ;   ii = ii + 1
75      WRITE(ldtxt(ii),*) 'usr_def_nam  : read the user defined namelist (namusr_def) in namelist_cfg'     ;   ii = ii + 1
76      WRITE(ldtxt(ii),*) '~~~~~~~~~~~ '                                                                   ;   ii = ii + 1
77      WRITE(ldtxt(ii),*) '   Namelist namusr_def : GYRE case'                                             ;   ii = ii + 1
78      WRITE(ldtxt(ii),*) '      GYRE used as Benchmark (=T)                      ln_bench  = ', ln_bench  ;   ii = ii + 1
79      WRITE(ldtxt(ii),*) '      inverse resolution & implied domain size         nn_GYRE   = ', nn_GYRE   ;   ii = ii + 1
80      WRITE(ldtxt(ii),*) '         jpiglo = 30*nn_GYRE+2                            jpiglo = ', kpi       ;   ii = ii + 1
81      WRITE(ldtxt(ii),*) '         jpjglo = 20*nn_GYRE+2                            jpjglo = ', kpj       ;   ii = ii + 1
82      WRITE(ldtxt(ii),*) '      number of model levels                           jpkglo    = ', kpk       ;   ii = ii + 1
83      !
84   END SUBROUTINE usr_def_nam
85
86
87   SUBROUTINE usr_def_hgr( plamt , plamu , plamv  , plamf  ,   &   ! geographic position (required)
88      &                    pphit , pphiu , pphiv  , pphif  ,   &   !
89      &                    kff   , pff_f , pff_t  ,            &   ! Coriolis parameter  (if domain not on the sphere)
90      &                    pe1t  , pe1u  , pe1v   , pe1f   ,   &   ! scale factors       (required)
91      &                    pe2t  , pe2u  , pe2v   , pe2f   ,   &   !
92      &                    ke1e2u_v      , pe1e2u , pe1e2v   )     ! u- & v-surfaces (if gridsize reduction is used in strait(s))
93      !!----------------------------------------------------------------------
94      !!                  ***  ROUTINE usr_def_hgr  ***
95      !!
96      !! ** Purpose :   user defined mesh and Coriolis parameter
97      !!
98      !! ** Method  :   set all intent(out) argument to a proper value
99      !!
100      !!                Here GYRE configuration :
101      !!          Rectangular mid-latitude domain
102      !!          - with axes rotated by 45 degrees
103      !!          - a constant horizontal resolution of 106 km
104      !!          - on a beta-plane
105      !!
106      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees)
107      !!              - define coriolis parameter at f-point if the domain in not on the sphere (on beta-plane)
108      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters)
109      !!              - define u- & v-surfaces (if gridsize reduction is used in some straits) (in m2)
110      !!----------------------------------------------------------------------
111      REAL(wp), DIMENSION(:,:), INTENT(out) ::   plamt, plamu, plamv, plamf   ! longitude outputs                     [degrees]
112      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pphit, pphiu, pphiv, pphif   ! latitude outputs                      [degrees]
113      INTEGER                 , INTENT(out) ::   kff                          ! =1 Coriolis parameter computed here, =0 otherwise
114      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pff_f, pff_t                 ! Coriolis factor at f-point                [1/s]
115      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1t, pe1u, pe1v, pe1f       ! i-scale factors                             [m]
116      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe2t, pe2u, pe2v, pe2f       ! j-scale factors                             [m]
117      INTEGER                 , INTENT(out) ::   ke1e2u_v                     ! =1 u- & v-surfaces computed here, =0 otherwise
118      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1e2u, pe1e2v               ! u- & v-surfaces (if reduction in strait)   [m2]
119      !
120      INTEGER  ::   ji, jj               ! dummy loop indices
121      REAL(wp) ::   zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1  , ze1deg, zf0 ! local scalars
122      REAL(wp) ::   zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme      !   -      -
123      !!-------------------------------------------------------------------------------
124      !
125      IF( nn_timing == 1 )  CALL timing_start('usr_def_hgr')
126      !
127      !     !==  beta-plane with regular grid-spacing and rotated domain ==!  (GYRE configuration)
128      !
129      IF(lwp) WRITE(numout,*)
130      IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
131      !
132      !                       !==  grid point position  ==!
133      !
134      zlam1 = -85._wp                           ! position of gridpoint (i,j) = (1,jpjglo)
135      zphi1 =  29._wp
136      !
137      ze1 = 106000._wp / REAL( nn_GYRE , wp )   ! gridspacing in meters
138      !
139      zsin_alpha = - SQRT( 2._wp ) * 0.5_wp     ! angle: 45 degrees
140      zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
141      ze1deg = ze1 / (ra * rad)
142      zlam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
143      zphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
144      !   
145      IF( ln_bench ) THEN     ! benchmark: forced the resolution to be 106 km
146         ze1 = 106000._wp     ! but keep (lat,lon) at the right nn_GYRE resolution
147         CALL ctl_warn( ' GYRE used as Benchmark: e1=e2=106km, no need to adjust rdt, ahm,aht ' )
148      ENDIF
149      IF( nprint==1 .AND. lwp )   THEN
150         WRITE(numout,*) 'ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
151         WRITE(numout,*) 'ze1deg', ze1deg, 'zlam0', zlam0, 'zphi0', zphi0
152      ENDIF
153      !   
154      DO jj = 1, jpj 
155         DO ji = 1, jpi 
156            zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5 
157            zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5 
158            !   
159            !glamt(i,j) longitude at T-point
160            !gphit(i,j) latitude at T-point 
161            plamt(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
162            pphit(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
163            !   
164            !glamu(i,j) longitude at U-point
165            !gphiu(i,j) latitude at U-point
166            plamu(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
167            pphiu(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
168            !   
169            !glamv(i,j) longitude at V-point
170            !gphiv(i,j) latitude at V-point
171            plamv(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
172            pphiv(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
173            !
174            !glamf(i,j) longitude at F-point
175            !gphif(i,j) latitude at F-point
176            plamf(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
177            pphif(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
178         END DO
179      END DO
180      !
181      !                       !== Horizontal scale factors ==! (in meters)
182      !                     
183      !                                         ! constant grid spacing
184      pe1t(:,:) =  ze1     ;      pe2t(:,:) = ze1
185      pe1u(:,:) =  ze1     ;      pe2u(:,:) = ze1
186      pe1v(:,:) =  ze1     ;      pe2v(:,:) = ze1
187      pe1f(:,:) =  ze1     ;      pe2f(:,:) = ze1
188      !
189      !                                         ! NO reduction of grid size in some straits
190      ke1e2u_v = 0                              !    ==>> u_ & v_surfaces will be computed in dom_ghr routine
191      !
192      !
193      !                       !==  Coriolis parameter  ==!
194      kff = 1                                            !  indicate not to compute ff afterward
195      !
196      zbeta = 2. * omega * COS( rad * zphi1 ) / ra       ! beta at latitude zphi1
197      !SF we overwrite zphi0 (south point in latitude) used just above to define pphif (value of zphi0=15.5190567531966)
198      !SF for computation of Coriolis we keep the parameter of Hazeleger, W., and S. S. Drijfhout, JPO 1998.
199      zphi0 = 15._wp                                     !  latitude of the most southern grid point 
200      zf0   = 2. * omega * SIN( rad * zphi0 )            !  compute f0 1st point south
201      !
202      pff_f(:,:) = ( zf0 + zbeta * ABS( pphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south)
203      pff_t(:,:) = ( zf0 + zbeta * ABS( pphit(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south)
204      !
205      IF(lwp) WRITE(numout,*) '                           beta-plane used. beta = ', zbeta, ' 1/(s.m)'
206      !
207      IF( nn_timing == 1 )  CALL timing_stop('usr_def_hgr')
208      !
209   END SUBROUTINE usr_def_hgr
210
211
212   SUBROUTINE usr_def_zgr()
213       ! subroutine for vertical grid
214   END SUBROUTINE usr_def_zgr
215
216 
217   SUBROUTINE usr_def_ini( pts )
218      !!----------------------------------------------------------------------
219      !!                   ***  ROUTINE usr_def_ini  ***
220      !!
221      !! ** Purpose :   Initialization of the dynamics and tracers
222      !!                Here GYRE configuration example : (double gyre with rotated domain)
223      !!
224      !! ** Method  : - set temprature field
225      !!              - set salinity   field
226      !!----------------------------------------------------------------------
227      USE dom_oce, ONLY : gdept_0, tmask
228      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pts   ! T & S data
229      !
230      INTEGER :: ji, jj, jk  ! dummy loop indices
231      !!----------------------------------------------------------------------
232      !
233      IF(lwp) WRITE(numout,*)
234      IF(lwp) WRITE(numout,*) 'usr_def_ini : analytical definition of initial state '
235      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   T and S profiles deduced from LEVITUS '
236      !
237      DO jk = 1, jpk
238         DO jj = 1, jpj
239            DO ji = 1, jpi
240               pts(ji,jj,jk,jp_tem) =  (  (  16. - 12. * TANH( (gdept_0(ji,jj,jk) - 400) / 700 ) )   &
241                    &           * (-TANH( (500-gdept_0(ji,jj,jk)) / 150. ) + 1.) / 2.               &
242                    &           + ( 15. * ( 1. - TANH( (gdept_0(ji,jj,jk)-50.) / 1500.) )        &
243                    &           - 1.4 * TANH((gdept_0(ji,jj,jk)-100.) / 100.)                    &
244                    &           + 7.  * (1500. - gdept_0(ji,jj,jk) ) / 1500.)                     &
245                    &           * (-TANH( (gdept_0(ji,jj,jk) - 500.) / 150.) + 1.) / 2.  ) * tmask(ji,jj,jk)
246
247               pts(ji,jj,jk,jp_sal) =  (  (  36.25 - 1.13 * TANH( (gdept_0(ji,jj,jk) - 305) / 460 ) )  &
248                    &         * (-TANH((500. - gdept_0(ji,jj,jk)) / 150.) + 1.) / 2              &
249                    &         + ( 35.55 + 1.25 * (5000. - gdept_0(ji,jj,jk)) / 5000.          &
250                    &         - 1.62 * TANH( (gdept_0(ji,jj,jk) - 60.  ) / 650. )             &
251                    &         + 0.2  * TANH( (gdept_0(ji,jj,jk) - 35.  ) / 100. )             &
252                    &         + 0.2  * TANH( (gdept_0(ji,jj,jk) - 1000.) / 5000.) )           &
253                    &         * (-TANH((gdept_0(ji,jj,jk) - 500.) / 150.) + 1.) / 2  ) * tmask(ji,jj,jk)
254            END DO
255         END DO
256      END DO
257      !   
258   END SUBROUTINE usr_def_ini
259
260   !!======================================================================
261END MODULE usrdef
Note: See TracBrowser for help on using the repository browser.