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.

domzgr_s.h90 in tags/nemo_dev_x3/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: tags/nemo_dev_x3/NEMO/OPA_SRC/DOM/domzgr_s.h90 @ 105

Last change on this file since 105 was 105, checked in by cvs2svn, 20 years ago
This commit was manufactured by cvs2svn to create tag 'nemo_dev_x3'.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 9.3 KB

Line
1	!!----------------------------------------------------------------------
2	!! * domzgr.s.h90 *
3	!!----------------------------------------------------------------------
4
5	#if defined key_s_coord
6	!!----------------------------------------------------------------------
7	!! 'key_s_coord' : s-coordinate
8	!!----------------------------------------------------------------------
9
10	SUBROUTINE zgr_s
11	!!----------------------------------------------------------------------
12	!! * ROUTINE zgr_s *
13	!!
14	!! ** Purpose : define the s-coordinate system
15	!!
16	!! ** Method : s-coordinate
17	!! The depth of model levels is defined as the product of an
18	!! analytical function by the local bathymetry, while the vertical
19	!! scale factors are defined as the product of the first derivative
20	!! of the analytical function by the bathymetry.
21	!! (this solution save memory as depth and scale factors are not
22	!! 3d fields)
23	!! - Read bathymetry (in meters) at t-point and compute the
24	!! bathymetry at u-, v-, and f-points.
25	!! hbatu = mi( hbatt )
26	!! hbatv = mj( hbatt )
27	!! hbatf = mi( mj( hbatt ) )
28	!! - Compute gsigt, gsigw, esigt, esigw from an analytical
29	!! function and its derivative given as statement function.
30	!! gsigt(k) = fssig (k+0.5)
31	!! gsigw(k) = fssig (k )
32	!! esigt(k) = fsdsig(k+0.5)
33	!! esigw(k) = fsdsig(k )
34	!! This routine is given as an example, it must be modified
35	!! following the user s desiderata. nevertheless, the output as
36	!! well as the way to compute the model levels and scale factors
37	!! must be respected in order to insure second order a!!uracy
38	!! schemes.
39	!!
40	!! Reference :
41	!! Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
42	!!
43	!! History :
44	!! ! 95-12 (G. Madec) Original code : s vertical coordinate
45	!! ! 97-07 (G. Madec) lbc_lnk call
46	!! ! 97-04 (J.-O. Beismann)
47	!! 8.5 ! 02-06 (G. Madec) F90: Free form and module
48	!!----------------------------------------------------------------------
49	!! * Local declarations
50	INTEGER :: in
51	INTEGER :: ji, jj, jk
52	REAL(wp) :: fssig, fsdsig, pfkk
53	REAL(wp) :: zh0, zdz1, zdzn, zd, za, zb, zc, zdepmi
54
55	!!----------------------------------------------------------------------
56	!! OPA 8.5, LODYC-IPSL (2002)
57	!!----------------------------------------------------------------------
58
59	! a. Sigma stretching coefficient
60	! hyperbolic stretching sigma= za(in-jk)^3+zb(in-jk)^2+zc(in-jk)+zd
61	! hyperbolic stretching sigma= za(jk)^3+zb(jk)^2+zc(jk)+zd
62	! calculate stretching coefficients (to ensure
63	! ss(sigma) =1 at sigma=1, and -1 at sigma=-1)
64	! sigma =1 at the surface (k=1) and =-1 at the bottom (k=jpk)
65	! the only free parameters are the depths of the first level
66	! dz1 and of the last level(zdzn) in meters (given as a positive
67	! quantity).
68
69	fssig(pfkk) =-0.5* ( ( ( za(pfkk-1.) + zb ) (pfkk-1.) + zc ) * (pfkk-1.) + zd -1. )
70
71	! b. Vertical derivative of sigma stretching coefficient
72
73	fsdsig(pfkk)= 0.5( ( 3.za(pfkk-1.) + 2.zb )* (pfkk-1.) + zc )
74
75	! ------------------ stretching polynomial ------------------------
76
77	! SPEM coefficients
78	! spem opa: kopa= jpk-kspem
79	! pfkk =in-jk (spem) <==> pfkk=in-jk-jpk=-jk-1 (opa)
80
81	zh0 = 5500.
82	zdz1 = 800.
83	zdzn = 17.5
84	in = 19
85
86	IF(lwp) THEN
87	WRITE(numout,*) ' stretching with a third order polynamial'
88	WRITE(numout,*) ' h0 = ',zh0
89	WRITE(numout,*) ' at ',zh0,' dz(1)= z(1)-z(0) = ',zdz1
90	WRITE(numout,*) ' dz(n)= z(n)-z(n-1) = ',zdzn
91	ENDIF
92
93	za = - 2.*(zdz1+zdzn)/zh0/float(in-1)/float(in-2) &
94	+ 4./float(in)/float(in-1)/float(in-2)
95
96	zb = - 6./float(in-1)/float(in-2) &
97	+ 2.(zdz1float(in+1)+ zdznfloat(2in-1)) &
98	/zh0/float(in-2)/float(in-1)
99
100	zc = -2.zdznfloat(in)/zh0/float(in-2) &
101	-2.zdz1float(in)/zh0/float(in-1)/float(in-2) &
102	+(6.* float(in) - 4.)/float(in)/float(in-1)/float(in-2)
103
104	zd = 1.
105
106	!!----------------------------------------------------------------------
107	!! OPA 8.5, LODYC-IPSL (2002)
108	!!----------------------------------------------------------------------
109
110
111	! 1. Lecture and computation of hbat fields
112	! -----------------------------------------
113
114	! 1.1 Read hbatt (in meters)
115	! 1.2 Set the hbatt to negative value ???
116
117	hbatt(:,:) = - bathy(:,:)
118
119
120	! 1.3 Set a minimum depth (zdepmi)
121
122	zdepmi = -50.
123	IF( zdepmi >= 0. ) THEN
124	IF(lwp) WRITE(numout,*)'domzgr: the minimum depth must be < 0'
125	STOP 'domzgr'
126	ENDIF
127	DO jj = 1, jpj
128	DO ji = 1, jpi
129	hbatt(ji,jj) = MIN( hbatt(ji,jj), zdepmi )
130	hbatu(ji,jj) = zdepmi
131	hbatv(ji,jj) = zdepmi
132	hbatf(ji,jj) = zdepmi
133	END DO
134	END DO
135
136	CALL lbc_lnk( hbatt, 'T', 1. )
137
138	IF(lwp) THEN
139	WRITE(numout,*)
140	WRITE(numout,*) ' domzgr: minimum depth set to : ',zdepmi
141	WRITE(numout,*)
142	ENDIF
143
144	! 1.4 Control print
145
146	IF(lwp) THEN
147	WRITE(numout,*)
148	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
149	WRITE(numout,*)
150	CALL prihre(hbatt(1,1),jpi,jpj,1,jpi,1,1,jpj,1,0.,numout)
151	ENDIF
152
153	! 1.5 Compute hbat fields at u-, v-, f-points
154
155	DO jj = 1, jpjm1
156	DO ji = 1, jpim1
157	hbatu(ji,jj)= 0.5 ( hbatt(ji ,jj)+hbatt(ji+1,jj ) ) umask(ji,jj,1)
158	hbatv(ji,jj)= 0.5 ( hbatt(ji ,jj)+hbatt(ji ,jj+1) ) vmask(ji,jj,1)
159	hbatf(ji,jj)= 0.25*( hbatt(ji ,jj)+hbatt(ji ,jj+1) &
160	+hbatt(ji+1,jj)+hbatt(ji+1,jj+1) ) * fmask(ji,jj,1)
161	END DO
162	END DO
163
164	CALL lbc_lnk( hbatu, 'U', 1. )
165	CALL lbc_lnk( hbatv, 'V', 1. )
166	CALL lbc_lnk( hbatf, 'F', 1. )
167
168
169	! 2. Computation of gsig and esig fields
170	! --------------------------------------
171
172	! 2.1 Coefficients for model level depth at w- and t-levels
173
174	DO jk = 1, jpk
175	gsigw(jk) = -fssig ( float(jk) )
176	gsigt(jk) = -fssig ( float(jk)+0.5)
177	END DO
178
179	! 2.2 Coefficients for vertical scale factors at w-, t- levels
180
181	DO jk = 1, jpk
182	esigw(jk)=fsdsig( float(jk) )
183	esigt(jk)=fsdsig( float(jk)+0.5)
184	END DO
185
186	! 2.3 Coefficients for vertical depth as the sum of e3w scale factors
187
188	gsi3w(1) = 0.5 * esigw(1)
189	DO jk = 2, jpk
190	gsi3w(jk) = gsi3w(jk-1)+ esigw(jk)
191	END DO
192
193	! 2.4 Control print
194
195	IF(lwp) THEN
196	WRITE(numout,*)
197	WRITE(numout,*) ' domzgr: vertical coefficients for model level'
198	WRITE(numout,9400)
199	WRITE(numout,9410) (jk,gsigt(jk),gsigw(jk),esigt(jk),esigw(jk),gsi3w(jk),jk=1,jpk)
200	ENDIF
201	9400 FORMAT(9x,' level gsigt gsigw esigt esigw gsi3w')
202	9410 FORMAT(10x,i4,5f9.2)
203
204	IF(lwp) THEN
205	WRITE(numout,*)
206	WRITE(numout,*) ' domzgr: vertical coordinates : point (10,10,k)'
207	WRITE(numout,*) ' ~~~~~~ --------------------'
208	WRITE(numout,9420)
209	WRITE(numout,9430) (jk,fsdept(10,10,jk),fsdepw(10,10,jk), &
210	fse3t (10,10,jk),fse3w (10,10,jk),jk=1,jpk)
211	ENDIF
212
213	9420 FORMAT(9x,' level gdept gdepw gde3w e3t e3w ')
214	9430 FORMAT(10x,i4,4f9.2)
215
216	DO jk = 1, jpk
217	DO jj = 1, jpj
218	DO ji = 1, jpi
219	IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
220	IF(lwp) THEN
221	WRITE(numout,*)
222	WRITE(numout,*) ' e r r o r : e3w or e3t =< 0 '
223	WRITE(numout,*) ' ========= --------------- '
224	WRITE(numout,*)
225	WRITE(numout,*) ' point (i,j,k)= ',ji,jj,jk
226	WRITE(numout,*)
227	ENDIF
228	STOP 'domzgr'
229	ENDIF
230	IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN
231	IF(lwp) THEN
232	WRITE(numout,*)
233	WRITE(numout,*) ' e r r o r : gdepw or gdept < 0 '
234	WRITE(numout,*) ' ========= ------------------ '
235	WRITE(numout,*)
236	WRITE(numout,*) ' point (i,j,k)= ',ji,jj,jk
237	WRITE(numout,*)
238	ENDIF
239	STOP 'domzgr'
240	ENDIF
241	END DO
242	END DO
243	END DO
244
245	END SUBROUTINE zgr_s
246
247	#else
248	!!----------------------------------------------------------------------
249	!! Default option : Empty routine
250	!!----------------------------------------------------------------------
251	SUBROUTINE zgr_s
252	END SUBROUTINE zgr_s
253	#endif

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format