1 | ! -*- Mode: f90 -*- |
---|
2 | PROGRAM cote |
---|
3 | USE declare |
---|
4 | USE modeles |
---|
5 | USE mod_lecdata |
---|
6 | USE formula |
---|
7 | USE math |
---|
8 | USE mod_rectifie |
---|
9 | USE mod_lbc |
---|
10 | USE mod_norma |
---|
11 | USE bords |
---|
12 | USE mod_inter |
---|
13 | USE mod_prih |
---|
14 | USE mod_wri_wei |
---|
15 | USE mod_inipar |
---|
16 | USE fliocom |
---|
17 | USE getincom |
---|
18 | USE errioipsl |
---|
19 | !!! |
---|
20 | IMPLICIT NONE |
---|
21 | !!! |
---|
22 | !!! A partir d'un fichier de poids, complete avec les rivieres et |
---|
23 | !!! les poids pour le run off |
---|
24 | !! On met a zero les poids sur les points non cotiers, et on renormalise |
---|
25 | !! |
---|
26 | INTEGER (kind=il), PARAMETER :: jpr = 52 ! Nombre de rivieres |
---|
27 | REAL (kind=rl) :: zsuma, zsumo |
---|
28 | REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: za |
---|
29 | REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: zo |
---|
30 | !! |
---|
31 | CHARACTER (LEN = 1) :: clriv |
---|
32 | CHARACTER (LEN = 30) :: clname |
---|
33 | INTEGER (kind=il) :: ja, jo, jo2, jo3, jn, jno, joi, joj, jai, jaj, jr, kn ! indices de boucle |
---|
34 | INTEGER (kind=il) :: jai_p, jai_m, jaj_p, jaj_m, n1 |
---|
35 | INTEGER (kind=il) :: joi_p, joi_m, joj_p, joj_m |
---|
36 | REAL (kind=rl) :: z_d_1, z_d_2, zzmin |
---|
37 | INTEGER (kind=il) :: nriv = 32, nlmd = 35 |
---|
38 | INTEGER (kind=il), DIMENSION (:, :), ALLOCATABLE :: jx ! Point trouve sous une maille atm |
---|
39 | INTEGER (kind=il) :: kom, ko |
---|
40 | LOGICAL, DIMENSION (:), ALLOCATABLE :: lacot, laland, laoce, lanoroute |
---|
41 | LOGICAL, DIMENSION (:), ALLOCATABLE :: locot, looce |
---|
42 | INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: ktemp |
---|
43 | REAL(kind=rl), DIMENSION (:), ALLOCATABLE :: wtemp, d_cote |
---|
44 | INTEGER (kind=il), DIMENSION (:,:), ALLOCATABLE :: ktemp2 |
---|
45 | REAL(kind=rl), DIMENSION (:,:), ALLOCATABLE :: wtemp2 |
---|
46 | INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: iamsk |
---|
47 | INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: iomsk |
---|
48 | INTEGER (kind=il) :: jovoid, jot, inum |
---|
49 | CHARACTER (len=80) :: cfname4, cfname8 |
---|
50 | LOGICAL :: ltrouve, ltrouve2 |
---|
51 | INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: itarget |
---|
52 | REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: wtarget |
---|
53 | INTEGER (kind=il), DIMENSION (:,:), ALLOCATABLE :: i2a |
---|
54 | INTEGER (kind=il) :: no, ntrouve |
---|
55 | INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: jo_n |
---|
56 | LOGICAL :: l_fulldiag = .TRUE. |
---|
57 | CHARACTER (len=180) :: c_date, c_time, c_zone, c_tmp |
---|
58 | CHARACTER (len=80) :: clarg |
---|
59 | INTEGER :: narg |
---|
60 | REAL(kind=rl), DIMENSION (:), ALLOCATABLE :: tmask_atl, tmask_noclose, tmask_nomed, tmask_med |
---|
61 | REAL(kind=rl), DIMENSION (:, :), ALLOCATABLE :: z2o |
---|
62 | CHARACTER (len=40) :: cl_bassin |
---|
63 | INTEGER (kind=il) :: il_ncid |
---|
64 | !! |
---|
65 | INTEGER (kind=il), DIMENSION (2) :: jind |
---|
66 | INTEGER (kind=il) :: ierr |
---|
67 | !! |
---|
68 | !! Recherche voisins proches |
---|
69 | INTEGER (kind=il) :: jp_nb |
---|
70 | INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: jo_tab |
---|
71 | INTEGER (kind=il) :: jo_nb |
---|
72 | !! |
---|
73 | INTEGER :: nid_it |
---|
74 | !! |
---|
75 | CHARACTER (len = 50) :: cldiag_a2o, clw_a2o, clw_a2o_cdf, clw_a2o_mct |
---|
76 | CHARACTER (len = 25) :: mo_name, a2o_name |
---|
77 | CHARACTER (len = 1) :: c_i, c_r |
---|
78 | |
---|
79 | !! Lecture de la ligne de comande |
---|
80 | narg = iargc() |
---|
81 | IF ( narg > 0 ) THEN |
---|
82 | CALL getarg ( 1, clarg) |
---|
83 | SELECT CASE (TRIM(clarg)) |
---|
84 | CASE Default |
---|
85 | IF ( INDEX ( TRIM(clarg), '.def' ) /= 0 ) THEN |
---|
86 | WRITE (nout,*) 'Lecture des parametres dans ', TRIM(clarg) |
---|
87 | CALL getin_name ( TRIM (clarg) ) |
---|
88 | ELSE |
---|
89 | WRITE (nout,*) 'Lecture des parametres dans run.def' |
---|
90 | ENDIF |
---|
91 | END SELECT |
---|
92 | END IF |
---|
93 | !! |
---|
94 | CALL inipar |
---|
95 | CALL alloc_modeles |
---|
96 | |
---|
97 | !! |
---|
98 | ALLOCATE (za (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'za', lreset=.TRUE., crout='cotes') |
---|
99 | ALLOCATE (zo (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'zo') |
---|
100 | ALLOCATE (jx (jpo2a,2) , STAT=ierr) ; CALL chk_allo (ierr, 'jx') |
---|
101 | ALLOCATE (lacot (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'lacot') ; lacot = .FALSE. |
---|
102 | ALLOCATE (laland (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'laland') ; laland = .FALSE. |
---|
103 | ALLOCATE (laoce (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'laoce') ; laoce = .FALSE. |
---|
104 | ! |
---|
105 | ALLOCATE (locot (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'locot') ; locot = .FALSE. |
---|
106 | ALLOCATE (looce (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'looce') ; looce = .FALSE. |
---|
107 | ! |
---|
108 | ALLOCATE (ktemp (jpa2o) , STAT=ierr) ; CALL chk_allo (ierr, 'ktemp') |
---|
109 | ALLOCATE (wtemp (jpa2o) , STAT=ierr) ; CALL chk_allo (ierr, 'wtemp') |
---|
110 | ALLOCATE (iamsk (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'iamsk') |
---|
111 | ALLOCATE (iomsk (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'iomsk') |
---|
112 | ALLOCATE (d_cote (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'd_cote') |
---|
113 | ! |
---|
114 | ALLOCATE (itarget (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'itarget') |
---|
115 | ALLOCATE (wtarget (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'wtarget') |
---|
116 | ALLOCATE (i2a (jpai, jpaj), STAT=ierr) ; CALL chk_allo (ierr, 'i2a') |
---|
117 | ! |
---|
118 | ALLOCATE (tmask_atl (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmask_atl') |
---|
119 | ALLOCATE (tmask_noclose (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmask_noclose') |
---|
120 | ALLOCATE (tmask_nomed (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmasl_nomed') |
---|
121 | ALLOCATE (tmask_med (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmask_med') |
---|
122 | ALLOCATE (z2o (jpoi, jpoj), STAT=ierr) ; CALL chk_allo (ierr, 'z2o') |
---|
123 | ! |
---|
124 | jp_nb = jpon |
---|
125 | ALLOCATE (jo_tab (jpon)) ; CALL chk_allo (ierr, 'jo_tab') |
---|
126 | ! |
---|
127 | IF (l_limit_iosize) l_fulldiag = .FALSE. |
---|
128 | !! |
---|
129 | !! Initialisation |
---|
130 | ngrd = 11 ; nsrf = 12 ; nmsk = 13 ; nchk = 14 ; ndeb = 9 ; nbug1 = 10 |
---|
131 | nwei4o2a = 81 ; nwei8o2a = 82 ; nwei4a2o = 83 ; nwei8a2o = 84 |
---|
132 | ! |
---|
133 | nout = 6 |
---|
134 | IF (nout /= 6 .AND. nout /= 0 ) & |
---|
135 | & OPEN (unit = nout, file = 'cote' // TRIM(c_suffix) // '.out', action = 'WRITE', status = 'unknown', position = 'REWIND' ) |
---|
136 | !! |
---|
137 | !! Diagnostics files |
---|
138 | !! |
---|
139 | IF (nout /= 6 ) & |
---|
140 | & OPEN (unit = nout, file = 'cotes' // TRIM(c_suffix) // '.out', position = 'REWIND', action = 'WRITE', status = 'replace' ) |
---|
141 | !! |
---|
142 | !! Read coordinates of all models, and weights |
---|
143 | !! |
---|
144 | CALL lecdata (lecweights = .TRUE., leco2amask = .TRUE. ) |
---|
145 | |
---|
146 | !! Lecture masques de bassin |
---|
147 | IF ( TRIM (cotyp) == 'orca4' .OR. & |
---|
148 | & TRIM (cotyp) == 'orca2' .OR. TRIM (cotyp) == 'orca2.0' .OR. TRIM (cotyp) == 'orca2.1' .OR. TRIM (cotyp) == 'orca2.3' .OR. & |
---|
149 | & TRIM (cotyp) == 'orca1' & |
---|
150 | & ) THEN |
---|
151 | cl_bassin = TRIM (cotyp) |
---|
152 | !-$$ IF (c_period /= 'none' ) cl_bassin = TRIM (cl_bassin) // TRIM (c_period) |
---|
153 | cl_bassin = TRIM (cl_bassin) // '.nc' |
---|
154 | WRITE (nout,*) 'Reading bassin masks in : ', TRIM (cl_bassin) |
---|
155 | CALL flioopfd (TRIM (cl_bassin), il_ncid ) |
---|
156 | CALL fliogetv (il_ncid, 'mask_atl' , z2o ) |
---|
157 | tmask_atl = RESHAPE (z2o, (/ jpon /) ) |
---|
158 | CALL fliogetv (il_ncid, 'mask_noclose', z2o ) |
---|
159 | tmask_noclose = RESHAPE (z2o, (/ jpon /) ) |
---|
160 | CALL fliogetv (il_ncid, 'mask_nomed' , z2o ) |
---|
161 | tmask_nomed = RESHAPE (z2o, (/ jpon /) ) |
---|
162 | CALL flioclo (il_ncid) |
---|
163 | ELSE |
---|
164 | tmask_atl = REAL (1-iomskt, rl) |
---|
165 | tmask_noclose = REAL (1-iomskt, rl) |
---|
166 | tmask_nomed = REAL (1-iomskt, rl) |
---|
167 | END IF |
---|
168 | ! |
---|
169 | WRITE (unit=nout, fmt=*) 'iomskt' |
---|
170 | CALL prihin (RESHAPE(iomskt ,(/jpoi,jpoj/)), 1_il, nout) |
---|
171 | WRITE (unit=nout, fmt=*) 'tmask_atl' |
---|
172 | CALL prihin (RESHAPE(NINT(tmask_atl) ,(/jpoi,jpoj/)), 1_il, nout) |
---|
173 | WRITE (unit=nout, fmt=*) 'tmask_med' |
---|
174 | !CALL prihin (RESHAPE(NINT(tmask_med) ,(/jpoi,jpoj/)), 1_il, nout) |
---|
175 | WRITE (unit=nout, fmt=*) 'tmask_noclose' |
---|
176 | !CALL prihin (RESHAPE(NINT(tmask_noclose),(/jpoi,jpoj/)), 1_il, nout) |
---|
177 | ! |
---|
178 | ! Set correct periodicity to grids |
---|
179 | CALL lbc (xolont, noperio, 'T', jpoi, jpoj) |
---|
180 | CALL lbc (xolonu, noperio, 'U', jpoi, jpoj) |
---|
181 | CALL lbc (xolonv, noperio, 'V', jpoi, jpoj) |
---|
182 | CALL lbc (xolonf, noperio, 'F', jpoi, jpoj) |
---|
183 | CALL lbc (xolatt, noperio, 'T', jpoi, jpoj) |
---|
184 | CALL lbc (xolatu, noperio, 'U', jpoi, jpoj) |
---|
185 | CALL lbc (xolatv, noperio, 'V', jpoi, jpoj) |
---|
186 | CALL lbc (xolatf, noperio, 'F', jpoi, jpoj) |
---|
187 | CALL lbc (xosrft, noperio, 'T', jpoi, jpoj) |
---|
188 | CALL lbc (iomskt, noperio, 'T', jpoi, jpoj) |
---|
189 | ! |
---|
190 | CALL lbc (xalont, naperio, 'T', jpai, jpaj) |
---|
191 | CALL lbc (xalonu, naperio, 'U', jpai, jpaj) |
---|
192 | CALL lbc (xalonv, naperio, 'V', jpai, jpaj) |
---|
193 | CALL lbc (xalatt, naperio, 'T', jpai, jpaj) |
---|
194 | CALL lbc (xalatu, naperio, 'U', jpai, jpaj) |
---|
195 | CALL lbc (xalatv, naperio, 'V', jpai, jpaj) |
---|
196 | CALL lbc (xasrft, naperio, 'T', jpai, jpaj) |
---|
197 | CALL lbc (iamskt, naperio, 'T', jpai, jpaj) |
---|
198 | !! |
---|
199 | CALL rectifie |
---|
200 | !! |
---|
201 | CALL fliocrfd ('itarget' // TRIM(c_suffix), (/'x', 'y'/), (/jpai, jpaj/), nid_it, mode='REPLACE') |
---|
202 | CALL flioputa (nid_it, "?", "title" , "itarget" ) |
---|
203 | CALL flioputa (nid_it, '?', 'Comment', TRIM(c_comment) ) |
---|
204 | CALL DATE_AND_TIME (c_date, c_time, c_zone ) |
---|
205 | CALL flioputa (nid_it, "?", "history" , "Created: "//c_date(1:4)//"-"//c_date(5:6)//"-"// & |
---|
206 | & c_date(7:8)//" "//c_time(1:2)//"h"//c_time(3:4)//" GMT"//TRIM(c_zone) ) |
---|
207 | CALL flioputa (nid_it, "?", "method" , "MOSAIC" ) |
---|
208 | CALL flioputa (nid_it, "?", "source_grid" , "curvilinear" ) |
---|
209 | CALL flioputa (nid_it, "?", "dest_grid" , "curvilinear" ) |
---|
210 | CALL flioputa (nid_it, "?", "Institution" , "IPSL" ) |
---|
211 | CALL flioputa (nid_it, "?", "Model" , "IPSL CM6" ) |
---|
212 | CALL GET_ENVIRONMENT_VARIABLE ( NAME="HOSTNAME", VALUE=c_tmp, TRIM_NAME=.TRUE., STATUS=ierr) |
---|
213 | IF ( ierr == 0 ) THEN |
---|
214 | CALL flioputa (nid_it, "?", "HOSTNAME" , TRIM(c_tmp) ) |
---|
215 | ELSE |
---|
216 | WRITE (nout,*) 'Environment variable not found : $HOSTNAME' |
---|
217 | END IF |
---|
218 | CALL GET_ENVIRONMENT_VARIABLE ( NAME="LOGNAME" , VALUE=c_tmp, TRIM_NAME=.TRUE., STATUS=ierr) |
---|
219 | IF ( ierr == 0 ) THEN |
---|
220 | CALL flioputa (nid_it, "?", "LOGNAME" , TRIM(c_tmp) ) |
---|
221 | ELSE |
---|
222 | WRITE (nout,*) 'Environment variable not found : $LOGNAME' |
---|
223 | END IF |
---|
224 | |
---|
225 | CALL fliodefv (nid_it, 'lon', (/1, 2/), units="degree_north", v_t=flio_r4 ) |
---|
226 | CALL fliodefv (nid_it, 'lat', (/1, 2/), units="degree_east", v_t=flio_r4 ) |
---|
227 | CALL fliodefv (nid_it, 'o2amask', (/1,2/), v_t=flio_r4) |
---|
228 | CALL fliodefv (nid_it, 'laland', (/1,2/), v_t=flio_i4) |
---|
229 | CALL fliodefv (nid_it, 'lacot' , (/1,2/), v_t=flio_i4) |
---|
230 | CALL fliodefv (nid_it, 'laoce' , (/1,2/), v_t=flio_i4) |
---|
231 | CALL fliodefv (nid_it, 'iamskp', (/1,2/), v_t=flio_i4) |
---|
232 | |
---|
233 | CALL flioputv (nid_it, 'lon' , RESHAPE ( xalont, (/jpai, jpaj/)) ) |
---|
234 | CALL flioputv (nid_it, 'lat' , RESHAPE ( xalatt, (/jpai, jpaj/)) ) |
---|
235 | CALL flioputv (nid_it, 'o2amask', RESHAPE (o2amask, (/jpai, jpaj/)) ) |
---|
236 | CALL flioputv (nid_it, 'iamskp' , RESHAPE (iamskp , (/jpai, jpaj/)) ) |
---|
237 | !! |
---|
238 | !! Compute edges of grid boxes |
---|
239 | !! |
---|
240 | CALL bord_a |
---|
241 | CALL bord_o |
---|
242 | CALL bord_oa |
---|
243 | !! |
---|
244 | mo_name = "mozaic" |
---|
245 | a2o_name = "a2o" |
---|
246 | !! |
---|
247 | WRITE (nout,*) 'o2amask (>0, <=0, >1) ', COUNT (o2amask > 0.0_rl), COUNT (o2amask <= 0.0_rl) & |
---|
248 | & , COUNT (o2amask >= 1.0_rl) |
---|
249 | i2a = 0_il |
---|
250 | WHERE (RESHAPE(o2amask,(/jpai,jpaj/)) .GT. 0.0_rl ) |
---|
251 | i2a = 1_il |
---|
252 | END WHERE |
---|
253 | WRITE (nout, *) 'i2a (o2amask > 0)' |
---|
254 | CALL prihin (i2a , 1_il, nout) |
---|
255 | !! |
---|
256 | !! Verif poids avant reduction |
---|
257 | !! |
---|
258 | WRITE (unit=nout,fmt=*) 'Verif poids globale' |
---|
259 | za = 0.0_rl |
---|
260 | WRITE (nout, *) 'Comptage xasrft ' , COUNT ( xasrft < 1.0E-10_rl ) |
---|
261 | WRITE (nout, *) 'Comptage o2amask ' , COUNT ( o2amask < 0.5_rl ) |
---|
262 | WHERE ( xasrft * o2amask /= 0.0_rl ) |
---|
263 | za = 1.0_rl / (xasrft * o2amask) |
---|
264 | END WHERE |
---|
265 | WRITE (nout, *) 'Somme za ', SUM ( za) |
---|
266 | WRITE (nout, *) 'Somme xasrft ', SUM ( xasrft) |
---|
267 | zsuma = DOT_PRODUCT (za, xasrft * o2amask) |
---|
268 | !! |
---|
269 | WRITE (nout, *) 'jma2o for coastal runoff : ', jma2o |
---|
270 | !! Atm mask interpolated to oce |
---|
271 | za = REAL (1_il-iamskt, KIND = rl ) |
---|
272 | CALL inter_a2o (za, a2omask ) |
---|
273 | !! 1 on Atm interpolated to oce |
---|
274 | za = REAL (1_il , KIND = rl ) |
---|
275 | CALL inter_a2o (za, a2ofull ) |
---|
276 | !! |
---|
277 | CALL verif ('Apres lecture') |
---|
278 | !! |
---|
279 | CALL bilan ("Poids ", c_case="complet", c_int_atm="Local", c_int_oce="Local") |
---|
280 | CALL bilan ("Poids ", c_case="ocean" , c_int_atm="Local", c_int_oce="Local") |
---|
281 | !! |
---|
282 | !! Remise des poids en integre |
---|
283 | DO jn = 1, jpa2o |
---|
284 | DO jo = 1, jpon |
---|
285 | ja = ka2o (jn, jo) |
---|
286 | IF (ja > 0) THEN |
---|
287 | IF ( o2amask (ja) .GT. 0.0_rl) & |
---|
288 | & wa2o (jn, jo) = wa2o (jn, jo) * xosrft (jo) / ( xasrft (ja) * o2amask (ja)) |
---|
289 | END IF |
---|
290 | END DO |
---|
291 | END DO |
---|
292 | ! |
---|
293 | CALL bilan ("Integre ", c_case="complet", c_int_atm="Integral", c_int_oce="Integral") |
---|
294 | CALL bilan ("Integre ", c_case="ocean" , c_int_atm="Integral", c_int_oce="Integral") |
---|
295 | CALL bilan ("Integre ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
296 | !! |
---|
297 | E_DRYRUN: IF (l_dryrun) THEN |
---|
298 | WRITE (nout,*) "Remplissage bidon des champs pour tester les I/O" |
---|
299 | jma2o = jma2or |
---|
300 | nvo = jma2or |
---|
301 | WRITE (nout,*) 'jma2o pour test : ', jma2o |
---|
302 | CALL RANDOM_SEED |
---|
303 | CALL RANDOM_NUMBER (wa2o) ; ka2o = NINT (wa2o*REAL(jpan,rl) ) |
---|
304 | |
---|
305 | DO jo = 1, jpon |
---|
306 | IF ( iomskt (jo) == 1 ) THEN |
---|
307 | wa2o (:, jo) = 0.0_rl |
---|
308 | ka2o (:, jo) = 0_il |
---|
309 | END IF |
---|
310 | END DO |
---|
311 | |
---|
312 | DO jo = 1, jpon |
---|
313 | DO jno = 1, jpa2o |
---|
314 | ja = ka2o (jno, jo) |
---|
315 | IF (ja > 0) THEN |
---|
316 | IF (iamskt (ja) == 1) THEN |
---|
317 | wa2o (jno, jo) = 0.0_rl ; ka2o (jno, jo) = 0_il |
---|
318 | END IF |
---|
319 | END IF |
---|
320 | END DO |
---|
321 | END DO |
---|
322 | ELSE |
---|
323 | L_lcoast: IF (lcoast) THEN |
---|
324 | !! Dataset 4 : idem 2, mais en interpolant seulement des points cotes atmosphere vers |
---|
325 | !! les points cotes ocean. |
---|
326 | WRITE (nout,*) 'Demarrage L_coast' |
---|
327 | !! |
---|
328 | WRITE (nout,*) 'Calculs points terre/oce/cote sur grille atmosphere' |
---|
329 | WHERE (o2amask .LT. epsfrac ) laland = .TRUE. |
---|
330 | WHERE (o2amask .GE. epsfrac .AND. o2amask .LE. (1.0_rl-epsfrac) ) lacot = .TRUE. |
---|
331 | WHERE ( o2amask .GT. (1.0_rl-epsfrac) ) laoce = .TRUE. |
---|
332 | |
---|
333 | WRITE (UNIT=nout, FMT=*) 'Points cotiers ', COUNT(lacot) |
---|
334 | |
---|
335 | WRITE (nout,*) 'Ecriture ' |
---|
336 | i2a = 0 |
---|
337 | WHERE (RESHAPE(laland,(/jpai,jpaj/))) i2a=1 |
---|
338 | CALL flioputv (nid_it, "laland", i2a) |
---|
339 | i2a = 0 |
---|
340 | WHERE (RESHAPE(laoce ,(/jpai,jpaj/))) i2a=1 |
---|
341 | CALL flioputv (nid_it, "laoce", i2a) |
---|
342 | i2a = 0 |
---|
343 | WHERE (RESHAPE(lacot ,(/jpai,jpaj/))) i2a=1 |
---|
344 | CALL flioputv (nid_it, "lacot", i2a) |
---|
345 | |
---|
346 | WRITE (nout,*) 'Determination des points cotes ocean ' |
---|
347 | ALLOCATE (ktemp2(jpoi,jpoj)) |
---|
348 | ktemp2 = 0_il |
---|
349 | locot = .FALSE. |
---|
350 | looce = .FALSE. |
---|
351 | IF ( lcoast ) THEN |
---|
352 | DO jo = 1, jpon |
---|
353 | joi = moi (jo) ; joj = moj (jo) ! Indices 2D |
---|
354 | joi_p = joi+1 ; joi_m = joi-1 |
---|
355 | joj_p = MIN(joj+1, jpoj) ; joj_m = MAX (joj-1, 1) |
---|
356 | IF ( iomskt (jo) == 0_il .AND. iomskp (jo) == 0_il ) THEN |
---|
357 | IF ( iomskt (m1o (joi_p, joj )) == 1 & |
---|
358 | .OR. iomskt (m1o (joi_m, joj )) == 1 & |
---|
359 | .OR. iomskt (m1o (joi , joj_p)) == 1 & |
---|
360 | .OR. iomskt (m1o (joi , joj_m)) == 1 & |
---|
361 | .OR. iomskt (m1o (joi_p, joj_p)) == 1 & |
---|
362 | .OR. iomskt (m1o (joi_p, joj_m)) == 1 & |
---|
363 | .OR. iomskt (m1o (joi_m, joj_p)) == 1 & |
---|
364 | .OR. iomskt (m1o (joi_m, joj_m)) == 1 ) THEN |
---|
365 | locot (jo) = .TRUE. |
---|
366 | ktemp2 (joi, joj) = 1_il |
---|
367 | ELSE |
---|
368 | looce (jo) = .TRUE. |
---|
369 | END IF |
---|
370 | END IF |
---|
371 | END DO |
---|
372 | ELSE |
---|
373 | locot (:) = ( iomskt == 0 ) |
---|
374 | ENDIF |
---|
375 | CALL lbc (locot, noperio, 'T', jpoi, jpoj) |
---|
376 | CALL lbc (looce, noperio, 'T', jpoi, jpoj) |
---|
377 | WRITE (nout,*) "iomskt" |
---|
378 | CALL prihin (RESHAPE(iomskt,(/jpoi,jpoj/)), 1, nout) |
---|
379 | WRITE (nout,*) "Locot" |
---|
380 | CALL prihin (ktemp2, 1, nout) |
---|
381 | DEALLOCATE (ktemp2) |
---|
382 | WRITE (nout,*) 'Verif coherence (tout doit etre nul)' |
---|
383 | WRITE (nout,*) COUNT ( laland .AND. lacot) |
---|
384 | WRITE (nout,*) COUNT ( lacot .AND. laoce) |
---|
385 | WRITE (nout,*) COUNT ( laoce .AND. laland) |
---|
386 | WRITE (nout,*) COUNT ( .NOT. ( laland .OR. lacot .OR. laoce) ) |
---|
387 | |
---|
388 | !-$$ iamsk = 0 ; WHERE (lacot) iamsk = 1 |
---|
389 | !-$$ WRITE(nout,*) 'Fraction lacot' |
---|
390 | !-$$ CALL prihin(RESHAPE(iamsk,(/jpai,jpaj/)),1_il,nout) |
---|
391 | !-$$ |
---|
392 | !-$$ iamsk = 0 ; WHERE (laland) iamsk = 1 |
---|
393 | !-$$ WRITE(nout,*) 'Points avec terre laland' |
---|
394 | !-$$ CALL prihin(RESHAPE(iamsk,(/jpai,jpaj/)),1_il,nout) |
---|
395 | !-$$ |
---|
396 | !-$$ iamsk = 0 ; WHERE (laoce) iamsk = 1 |
---|
397 | !-$$ WRITE(nout,*) 'Fraction laoce' |
---|
398 | !-$$ CALL prihin(RESHAPE(iamsk,(/jpai,jpaj/)),1_il,nout) |
---|
399 | !-$$ |
---|
400 | !-$$ WRITE(nout,*) 'Periodicite atm' |
---|
401 | !-$$ CALL prihin(RESHAPE(iamskp,(/jpai,jpaj/)),1_il,nout) |
---|
402 | !! |
---|
403 | WRITE (nout,*) "iamskt" |
---|
404 | CALL prihin (RESHAPE(iamskt,(/jpai,jpaj/)), 1, nout) |
---|
405 | WRITE (nout,*) "Lacot" |
---|
406 | ALLOCATE (ktemp2 (jpai, jpaj)) |
---|
407 | ktemp2 = 0_il ; WHERE (RESHAPE(lacot, (/jpai, jpaj/))) ktemp2 = 1_il |
---|
408 | CALL prihin (ktemp2, 1, nout) |
---|
409 | !! |
---|
410 | WRITE (unit = nout, fmt = *) 'Nombre de points cotes atmosphere : ', COUNT (lacot) |
---|
411 | WRITE (unit = nout, fmt = *) 'Nombre de points cotes ocean : ', COUNT (locot) |
---|
412 | !! |
---|
413 | CALL bilan ("Avec cote ", c_case="complet", c_int_atm="Integral", c_int_oce="Integral") |
---|
414 | CALL bilan ("Avec cote ", c_case="ocean" , c_int_atm="Integral", c_int_oce="Integral") |
---|
415 | CALL bilan ("Avec cote ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
416 | ! |
---|
417 | ! Mise a zero des poids sur les points non cotiers |
---|
418 | DO jo = 1, jpon |
---|
419 | IF (.NOT. locot (jo) ) THEN |
---|
420 | wa2o (:, jo) = 0.0e0_rl |
---|
421 | ka2o (:, jo) = 0_il |
---|
422 | ELSE |
---|
423 | DO jn = 1, jpa2o |
---|
424 | ja = ka2o (jn, jo) |
---|
425 | IF (ja > 0 ) THEN |
---|
426 | IF (.NOT. lacot (ja) ) THEN |
---|
427 | wa2o (jn, jo) = 0.0e0_rl |
---|
428 | ka2o (jn, jo) = 0_il |
---|
429 | END IF |
---|
430 | END IF |
---|
431 | END DO |
---|
432 | END IF |
---|
433 | END DO |
---|
434 | CALL verif ('Apres remise a zero') |
---|
435 | !! |
---|
436 | !! Elimination des points redondant par periodicite |
---|
437 | DO jo = 1, jpon |
---|
438 | IF (iomskp (jo) .EQ. 1_il ) THEN |
---|
439 | wa2o (:,jo) = 0.0e0_rl |
---|
440 | ka2o (:,jo) = 0_il |
---|
441 | END IF |
---|
442 | END DO |
---|
443 | !! |
---|
444 | !! Controle |
---|
445 | !! |
---|
446 | CALL verif ('Apres elimination de periodicite') |
---|
447 | !! |
---|
448 | CALL bilan ("Apres remise a zero ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
449 | !! |
---|
450 | !! Renormalize |
---|
451 | ! Cas 1 |
---|
452 | RenormA : DO ja = 1, jpan |
---|
453 | IF (.NOT. lacot (ja) ) CYCLE RenormA |
---|
454 | z_d_1 = 0.0e0_rl |
---|
455 | kom = 0 |
---|
456 | RenormO : DO jo = 1, jpon |
---|
457 | IF (.NOT. locot (jo) ) CYCLE RenormO |
---|
458 | DO jn = 1, jpa2o |
---|
459 | IF ( ka2o (jn,jo) == ja ) THEN |
---|
460 | z_d_1 = z_d_1 + wa2o (jn,jo) |
---|
461 | kom = kom + 1 |
---|
462 | jx (kom,:) = (/ jn, jo /) |
---|
463 | !!$ WRITE (unit = nout, fmt = '(6I6,16F18.6) ') ja, jo, jn, kom, jx (kom, :), & |
---|
464 | !!$ wa2o (jn,jo), wa2o (jx (kom, 1), jx (kom, 2)), z_d_1, & |
---|
465 | !!$ xosrft (jo), xasrft (ja) |
---|
466 | END IF |
---|
467 | END DO |
---|
468 | ENDDO RenormO |
---|
469 | IF (kom /= 0 ) THEN |
---|
470 | IF (z_d_1 == 0.0_rl ) THEN |
---|
471 | WRITE (nout, *) 'z_d_1 nul ', ja, xasrft (ja), lacot (ja), iamskt (ja), & |
---|
472 | & wa2o (:,ja), ka2o (:,ja) |
---|
473 | STOP |
---|
474 | ELSE |
---|
475 | DO ko = 1, kom |
---|
476 | wa2o (jx (ko,1), jx (ko,2)) = wa2o (jx (ko,1), jx (ko,2)) / z_d_1 |
---|
477 | END DO |
---|
478 | ENDIF |
---|
479 | END IF |
---|
480 | ENDDO RenormA |
---|
481 | !! |
---|
482 | !! Calcul des distances à la cote des points oceans |
---|
483 | !-$$ WRITE (nout,*) 'Distance cotes' |
---|
484 | !-$$ d_cote = rpi |
---|
485 | !-$$ Cal_dist_1: DO jo = 1, jpon |
---|
486 | !-$$ IF ( iomskt (jo) == 1 .OR. locot (jo) ) THEN |
---|
487 | !-$$ d_cote (jo) = 0.0E0_rl |
---|
488 | !-$$ ELSE |
---|
489 | !-$$ Cal_dist_2:DO jo2 = 1, jpon |
---|
490 | !-$$ IF ( locot (jo2) ) THEN |
---|
491 | !-$$ d_cote (jo) = MIN (d_cote (jo), REAL (geodist(xolont(jo), xolatt(jo), xolont(jo2), xolatt(jo2)), KIND=rl)) |
---|
492 | !-$$ END IF |
---|
493 | !-$$ END DO Cal_dist_2 |
---|
494 | !-$$ END IF |
---|
495 | !-$$ END DO Cal_dist_1 |
---|
496 | !-$$ d_cote = d_cote * r_earth |
---|
497 | !CALL prihre (RESHAPE(d_cote,(/jpoi,jpoj/)), 1.0E-5_rl, nout) |
---|
498 | WRITE (nout,*) 'Distance maxi a la cote (ocean) : ', MAXVAL (d_cote) |
---|
499 | CALL fliocrfd ("dist_cote" // TRIM(c_suffix), (/ "x", "y"/), (/jpoi, jpoj/), n1, MODE="REPLACE") |
---|
500 | CALL flioputa (n1, '?', 'Comment', TRIM(c_comment) ) |
---|
501 | CALL fliodefv (n1, "dist_cote", (/1, 2/)) |
---|
502 | CALL flioputv (n1, "dist_cote", RESHAPE (d_cote, (/jpoi, jpoj/) )) |
---|
503 | CALL flioclo (n1) |
---|
504 | ! |
---|
505 | !! Controle |
---|
506 | !! |
---|
507 | CALL verif ('Premiere normalisation') |
---|
508 | CALL gather_wei ('Premiere etape') |
---|
509 | CALL verif ('Premier gather') |
---|
510 | !! |
---|
511 | CALL bilan ("Renorm ", c_case="complet", c_int_atm="Integral", c_int_oce="Integral") |
---|
512 | CALL bilan ("Renorm ", c_case="ocean" , c_int_atm="Integral", c_int_oce="Integral") |
---|
513 | CALL bilan ("Renorm ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
514 | !! |
---|
515 | !! Traitement des points atmospheres non atteint |
---|
516 | !! |
---|
517 | total: IF (ltotal) THEN |
---|
518 | WRITE (nout,*) 'Va chercher tout les points terre (ltotal=TRUE)' |
---|
519 | SeekAtm_1_1: DO ja = 1, jpan |
---|
520 | !-$$ WRITE(nout,*) 'Point atm : ', ja, m2ai(ja), m2aj(ja), itarget(ja), iamskp(ja), laland(ja), lacot(ja), laoce(ja) |
---|
521 | IF (itarget (ja) .GE. 1_il) CYCLE SeekAtm_1_1 |
---|
522 | IF (iamskp (ja) .EQ. 1_il) CYCLE SeekAtm_1_1 |
---|
523 | IF (.NOT. laland (ja) ) CYCLE SeekAtm_1_1 |
---|
524 | zzmin = REAL (cartdist(0.0_rl, 90.0_rl, 0.0_rl, -90.0_rl), KIND=rl) |
---|
525 | ltrouve = .FALSE. ; jot = 0_il |
---|
526 | !-$$ WRITE(nout,*) 'Point atm : ', ja, m2ai(ja), m2aj(ja) |
---|
527 | SeekOce_1_1: DO jo = 1, jpon ! Cherche point ocean le plus proche |
---|
528 | IF (.NOT. locot (jo) ) CYCLE SeekOce_1_1 |
---|
529 | IF (iomskp(jo) .EQ. 1_il) CYCLE SeekOce_1_1 |
---|
530 | z_d_1 = REAL (cartdist(xalont(ja), xalatt(ja), xolont(jo), xolatt(jo)), KIND=rl) |
---|
531 | IF (z_d_1 .LT. zzmin) THEN |
---|
532 | zzmin = z_d_1 |
---|
533 | ltrouve = .TRUE. ; jot = jo |
---|
534 | END IF |
---|
535 | ENDDO SeekOce_1_1 |
---|
536 | Trouve_1_1: IF (ltrouve) THEN |
---|
537 | ! Calcul du poids |
---|
538 | z_d_1 = 1.0_rl |
---|
539 | ! Rangement des poids |
---|
540 | IF (nvo (jot) < jpa2o ) THEN ! Reste un poids libre |
---|
541 | nvo (jot) = nvo (jot) + 1 |
---|
542 | wa2o (nvo (jot), jot) = z_d_1 |
---|
543 | ka2o (nvo (jot), jot) = ja |
---|
544 | !-$$ WRITE (nout,*) ' Vers point oce ', jot, m2oi(jot), m2oj(jot), nvo(jot) |
---|
545 | ELSE ! Trouve12 |
---|
546 | WRITE(nout,*) ' Plus de place ', ja, jot, nvo (jot), z_d_1 |
---|
547 | STOP |
---|
548 | ENDIF |
---|
549 | ELSE ! Trouve |
---|
550 | !! Message d'erreur |
---|
551 | WRITE (unit = nout, fmt = '("Pas de voisin proche pour ", 1I5,2I4)' ) ja, m2ai(ja), m2aj(ja) |
---|
552 | END IF Trouve_1_1 |
---|
553 | END DO SeekAtm_1_1 |
---|
554 | !! |
---|
555 | CALL verif ('ltotal avant gather') |
---|
556 | CALL bilan ("Ltotal avant gather ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral") |
---|
557 | CALL bilan ("Ltotal avant gather ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral") |
---|
558 | CALL bilan ("Ltotal avant gather ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
559 | !! |
---|
560 | CALL gather_wei ('apres total') |
---|
561 | CALL verif ( 'Cas total', l_non_a=.TRUE.) |
---|
562 | !! |
---|
563 | WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o |
---|
564 | !! |
---|
565 | CALL bilan ("Fin ltotal ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral") |
---|
566 | CALL bilan ("Fin ltotal ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral") |
---|
567 | CALL bilan ("Fin ltotal ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
568 | !! |
---|
569 | END IF Total |
---|
570 | !! |
---|
571 | Total_dist: IF (ltotal_dist) THEN |
---|
572 | !! On limite la recherche à une certaine distance de la cote. |
---|
573 | WRITE (nout,*) 'Cas ltotal_dist' |
---|
574 | SeekAtm_2_1: DO ja = 1, jpan |
---|
575 | IF (itarget (ja) .GE. 1_il) CYCLE SeekAtm_2_1 |
---|
576 | IF (iamskp (ja) .EQ. 1_il) CYCLE SeekAtm_2_1 |
---|
577 | IF (.NOT. laland (ja) ) CYCLE SeekAtm_2_1 |
---|
578 | zzmin = dist_max |
---|
579 | ltrouve = .FALSE. ; jot = 0_il |
---|
580 | !$$$ WRITE(nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja) |
---|
581 | SeekOce_2_1: DO jo = 1, jpon ! Cherche point ocean le plus proche |
---|
582 | IF (.NOT. locot (jo) ) CYCLE SeekOce_2_1 |
---|
583 | IF (iomskp(jo) .EQ. 1_il) CYCLE SeekOce_2_1 |
---|
584 | z_d_1 = r_earth * REAL (geodist(xalont(ja), xalatt(ja), xolont(jo), xolatt(jo)), KIND=rl) |
---|
585 | IF (z_d_1 .LT. zzmin) THEN |
---|
586 | zzmin = z_d_1 |
---|
587 | ltrouve = .TRUE. ; jot = jo |
---|
588 | END IF |
---|
589 | ENDDO SeekOce_2_1 |
---|
590 | Trouve_2_1: IF (ltrouve) THEN |
---|
591 | ! Calcul du poids |
---|
592 | z_d_1 = 1.0_rl |
---|
593 | IF (nvo (jot) < jpa2o ) THEN ! Reste un poids libre |
---|
594 | nvo (jot) = nvo (jot) + 1 |
---|
595 | wa2o (nvo (jot), jot) = z_d_1 |
---|
596 | ka2o (nvo (jot), jot) = ja |
---|
597 | !$$$ WRITE (nout,*) ' Vers point oce ', jot, m2oi(jot), m2oj(jot), nvo(jot) |
---|
598 | ELSE |
---|
599 | WRITE(nout,*) ' Plus de place ', ja, jot, nvo (jot), z_d_1 |
---|
600 | STOP |
---|
601 | ENDIF |
---|
602 | END IF Trouve_2_1 |
---|
603 | END DO SeekAtm_2_1 |
---|
604 | !! |
---|
605 | CALL verif ('Ltotal_dist avant gather') |
---|
606 | CALL bilan ("Ltotal_dist avant gather ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral") |
---|
607 | CALL bilan ("Ltotal_dist avant gather ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral") |
---|
608 | CALL bilan ("Ltotal_dist avant gather ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
609 | !! |
---|
610 | CALL gather_wei ('apres total_dist') |
---|
611 | CALL verif ( 'Cas total_dist') |
---|
612 | !! |
---|
613 | WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o |
---|
614 | !! |
---|
615 | CALL bilan ("Fin ltotal_dist ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral") |
---|
616 | CALL bilan ("Fin ltotal_dist ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral") |
---|
617 | CALL bilan ("Fin ltotal_dist ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
618 | !! |
---|
619 | END IF Total_dist |
---|
620 | !! |
---|
621 | Total_dist_2: IF (ltotal_dist_2) THEN |
---|
622 | !! On limite la recherche à une certaine distance de la cote. |
---|
623 | !! On etale sur les point oceans cotiers proches |
---|
624 | WRITE (nout,*) '-------------------' |
---|
625 | WRITE (nout,*) 'Cas ltotal_dist_2' |
---|
626 | WRITE (nout,*) 'Distance recherche : ', dist_max/1.0E3_rl, 'km' |
---|
627 | WRITE (nout,*) 'Distance etalement ocean : ', dist_max_voisin/1.0E3_rl, 'km' |
---|
628 | !! |
---|
629 | SeekAtm_3_1: DO ja = 1, jpan |
---|
630 | IF (itarget (ja) .GE. 1_il) CYCLE SeekAtm_3_1 |
---|
631 | IF (iamskp (ja) .EQ. 1_il) CYCLE SeekAtm_3_1 |
---|
632 | IF (.NOT. laland (ja) ) CYCLE SeekAtm_3_1 |
---|
633 | zzmin = dist_max |
---|
634 | ltrouve = .FALSE. ; jot = 0_il ; ntrouve = 0 |
---|
635 | !WRITE(nout,*) 'Point atm : ', ja, m2ai(ja), m2aj(ja) |
---|
636 | SeekOce_3_1: DO jo = 1, jpon ! Cherche point ocean le plus proche |
---|
637 | IF (.NOT. locot (jo) ) CYCLE SeekOce_3_1 |
---|
638 | IF (iomskp (jo) .EQ. 1_il) CYCLE SeekOce_3_1 |
---|
639 | z_d_1 = r_earth * REAL (geodist (xalont (ja), xalatt (ja), xolont (jo), xolatt (jo)), KIND=rl) |
---|
640 | IF (z_d_1 .LT. zzmin) THEN |
---|
641 | zzmin = z_d_1 ; ltrouve = .TRUE. ; jot = jo ; ntrouve = ntrouve+1 |
---|
642 | END IF |
---|
643 | ENDDO SeekOce_3_1 |
---|
644 | WRITE(nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja), jot, ntrouve |
---|
645 | Trouve_3_1: IF (ltrouve) THEN |
---|
646 | ! On cherche les points proches du point ocean trouve |
---|
647 | z_d_1 = 0.0_rl ; jo_tab = 0 ; jo_nb = 0 |
---|
648 | SeekOce_3_2: DO jo2 = 1, jpon |
---|
649 | IF (.NOT. locot (jo2) ) CYCLE SeekOce_3_2 |
---|
650 | IF (iomskp (jo2) .EQ. 1_il) CYCLE SeekOce_3_2 |
---|
651 | IF ( r_earth * REAL (geodist (xolont (jot), xolatt (jot), xolont (jo2), xolatt (jo2)), KIND=rl) & |
---|
652 | & < dist_max_voisin) THEN |
---|
653 | IF (jo_nb >= jp_nb) THEN |
---|
654 | WRITE (nout,*) 'Trop de point proches ', ja, jot, jo_tab |
---|
655 | STOP |
---|
656 | ELSE |
---|
657 | jo_nb = jo_nb + 1 |
---|
658 | z_d_1 = z_d_1 + xosrft (jo2) |
---|
659 | jo_tab (jo_nb) = jo2 |
---|
660 | END IF |
---|
661 | END IF |
---|
662 | ! Calcul du poids |
---|
663 | END DO SeekOce_3_2 |
---|
664 | IF (jo_nb == 0) THEN |
---|
665 | WRITE (nout,*) 'Pas de voisin au voisin ! ', ja, jot, m2oi(jot), m2oj(jot), xolont(jot), xolatt(jot) |
---|
666 | WRITE (nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja) |
---|
667 | WRITE (nout,*) r_earth * REAL (geodist (xolont (jot), xolatt (jot), xolont (jot), xolatt (jot)), KIND=rl) |
---|
668 | STOP |
---|
669 | ELSE |
---|
670 | !-$$ WRITE(nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja), jot, jo_nb, jo_tab (1:jo_nb) |
---|
671 | DO jn = 1, jo_nb |
---|
672 | jo3 = jo_tab (jn) |
---|
673 | IF (nvo (jo3) < jpa2o ) THEN ! Reste un poids libre |
---|
674 | nvo (jo3) = nvo (jo3) + 1 |
---|
675 | wa2o (nvo (jo3), jo3) = xosrft (jo3) / z_d_1 |
---|
676 | ka2o (nvo (jo3), jo3) = ja |
---|
677 | !$$$ WRITE (nout,*) ' Vers point oce ', jot, m2oi(jot), m2oj(jot), nvo(jot) |
---|
678 | ELSE |
---|
679 | WRITE(nout,*) ' Plus de place ', ja, jo3, nvo (jo3), z_d_1 |
---|
680 | STOP |
---|
681 | ENDIF |
---|
682 | END DO |
---|
683 | END IF |
---|
684 | !-$$ ELSE |
---|
685 | !-$$ WRITE (nout,*) 'Pas de voisins proche trouve ' |
---|
686 | !-$$ WRITE (nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja) |
---|
687 | !-$$ WRITE (nout,*) dist_max, zzmin, z_d_1 |
---|
688 | !-$$ !STOP |
---|
689 | END IF Trouve_3_1 |
---|
690 | END DO SeekAtm_3_1 |
---|
691 | !! |
---|
692 | CALL verif ('Ltotal_dist_2 avant gather') |
---|
693 | CALL bilan ("Ltotal_dist_2 avant gather ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral") |
---|
694 | CALL bilan ("Ltotal_dist_2 avant gather ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral") |
---|
695 | CALL bilan ("Ltotal_dist_2 avant gather ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
696 | !! |
---|
697 | CALL gather_wei ('apres total_dist_2') |
---|
698 | CALL verif ( 'Cas total_dist_2') |
---|
699 | !! |
---|
700 | WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o |
---|
701 | !! |
---|
702 | CALL bilan ("Fin ltotal_dist_2 ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral") |
---|
703 | CALL bilan ("Fin ltotal_dist_2 ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral") |
---|
704 | CALL bilan ("Fin ltotal_dist_2 ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
705 | !! |
---|
706 | END IF Total_dist_2 |
---|
707 | !! |
---|
708 | Total_dist_3: IF (ltotal_dist_3) THEN |
---|
709 | !! On limite la recherche à une certaine distance de la cote. |
---|
710 | !! On etale sur les point oceans proches, cote ou pas |
---|
711 | WRITE (nout,*) '-------------------' |
---|
712 | WRITE (nout,*) 'Cas ltotal_dist_3' |
---|
713 | WRITE (nout,*) 'Distance recherche : ', dist_max/1.0E3_rl, 'km' |
---|
714 | WRITE (nout,*) 'Distance etalement ocean : ', dist_max_voisin/1.0E3_rl, 'km' |
---|
715 | !! |
---|
716 | SeekAtm_4_1: DO ja = 1, jpan |
---|
717 | IF (itarget (ja) .GE. 1_il) CYCLE SeekAtm_4_1 |
---|
718 | IF (iamskp (ja) .EQ. 1_il) CYCLE SeekAtm_4_1 |
---|
719 | IF (.NOT. laland (ja) ) CYCLE SeekAtm_4_1 |
---|
720 | zzmin = dist_max |
---|
721 | ltrouve = .FALSE. ; jot = 0_il |
---|
722 | !-$$ WRITE(nout,*) 'Point atm : ', ja, m2ai(ja), m2aj(ja) |
---|
723 | SeekOce_4_1: DO jo = 1, jpon ! Cherche point ocean le plus proche |
---|
724 | IF (iomskp (jo) .EQ. 1_il .OR. iomskt (jo) .EQ. 1_il ) CYCLE SeekOce_4_1 |
---|
725 | z_d_1 = r_earth * REAL (geodist (xalont (ja), xalatt (ja), xolont (jo), xolatt (jo)), KIND=rl) |
---|
726 | IF (z_d_1 .LT. zzmin) THEN |
---|
727 | zzmin = z_d_1 ; ltrouve = .TRUE. ; jot = jo |
---|
728 | END IF |
---|
729 | ENDDO SeekOce_4_1 |
---|
730 | !-$$ WRITE(nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja), jot |
---|
731 | Trouve_4_1: IF (ltrouve) THEN |
---|
732 | ! On cherche les points proches du point ocean trouve, le long de la cote et vers le large |
---|
733 | z_d_1 = 0.0_rl ; jo_tab = 0 ; jo_nb = 0 |
---|
734 | SeekOce_4_2: DO jo2 = 1, jpon |
---|
735 | IF (iomskp (jo2) .EQ. 1_il .OR. iomskt (jo2) .EQ. 1_il) CYCLE SeekOce_4_2 |
---|
736 | IF ( .NOT. l_large .AND. .NOT. locot (jo2)) CYCLE SeekOce_4_2 |
---|
737 | z_d_2 = r_earth * REAL (geodist (xolont (jot), xolatt (jot), xolont (jo2), xolatt (jo2)), KIND=rl) |
---|
738 | IF ( ( locot (jo2) .AND. z_d_2 < dist_max_cote ) & |
---|
739 | .OR. ( .NOT. locot (jo2) .AND. z_d_2 < dist_max_large ) ) THEN |
---|
740 | IF (jo_nb > jp_nb) THEN |
---|
741 | WRITE (nout,*) 'Trop de point proches ', ja, jot, jo_tab |
---|
742 | STOP |
---|
743 | ELSE |
---|
744 | jo_nb = jo_nb + 1 |
---|
745 | z_d_1 = z_d_1 + xosrft (jo2) |
---|
746 | jo_tab (jo_nb) = jo2 |
---|
747 | END IF |
---|
748 | END IF |
---|
749 | ! Calcul du poids |
---|
750 | END DO SeekOce_4_2 |
---|
751 | IF (jo_nb == 0) THEN |
---|
752 | WRITE (nout,*) 'Pas de voisin au voisin ! ', ja, jot, m2oi(jot), m2oj(jot), xolont(jot), xolatt(jot) |
---|
753 | WRITE (nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja) |
---|
754 | WRITE (nout,*) r_earth * REAL (geodist (xolont (jot), xolatt (jot), xolont (jot), xolatt (jot)), KIND=rl) |
---|
755 | STOP |
---|
756 | ELSE |
---|
757 | !-$$ WRITE(nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja), jot, jo_nb, jo_tab (1:jo_nb) |
---|
758 | DO jn = 1, jo_nb |
---|
759 | jo3 = jo_tab (jn) |
---|
760 | IF (nvo (jo3) < jpa2o ) THEN ! Reste un poids libre |
---|
761 | nvo (jo3) = nvo (jo3) + 1 |
---|
762 | wa2o (nvo (jo3), jo3) = xosrft (jo3) / z_d_1 |
---|
763 | ka2o (nvo (jo3), jo3) = ja |
---|
764 | !$$$ WRITE (nout,*) ' Vers point oce ', jot, m2oi(jot), m2oj(jot), nvo(jot) |
---|
765 | ELSE |
---|
766 | WRITE(nout,*) ' Plus de place ', ja, jo3, nvo (jo3), z_d_1 |
---|
767 | STOP |
---|
768 | ENDIF |
---|
769 | END DO |
---|
770 | END IF |
---|
771 | ELSE |
---|
772 | WRITE (nout,*) 'Pas de voisins proche trouve ' |
---|
773 | WRITE (nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja) |
---|
774 | WRITE (nout,*) dist_max, zzmin, z_d_1 |
---|
775 | STOP |
---|
776 | END IF Trouve_4_1 |
---|
777 | END DO SeekAtm_4_1 |
---|
778 | !! |
---|
779 | CALL verif ('Ltotal_dist_3 avant gather') |
---|
780 | CALL bilan ("Ltotal_dist_3 avant gather ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral") |
---|
781 | CALL bilan ("Ltotal_dist_3 avant gather ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral") |
---|
782 | CALL bilan ("Ltotal_dist_3 avant gather ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
783 | !! |
---|
784 | CALL gather_wei ('apres total_dist_3') |
---|
785 | CALL verif ( 'Cas total_dist_3') |
---|
786 | !! |
---|
787 | WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o |
---|
788 | !! |
---|
789 | CALL bilan ("Fin ltotal_dist_3 ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral") |
---|
790 | CALL bilan ("Fin ltotal_dist_3 ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral") |
---|
791 | CALL bilan ("Fin ltotal_dist_3 ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
792 | !! |
---|
793 | END IF Total_dist_3 |
---|
794 | !! |
---|
795 | Essai: IF (lessai) THEN |
---|
796 | WRITE (nout,*) '--- Methode essai ----' |
---|
797 | no = COUNT ( laland .AND. (itarget == 0_il)) |
---|
798 | WRITE (nout,*) 'Nombre de points atm non vises : ', no |
---|
799 | WRITE (nout,*) 'jma2o, jpa2o, jma2o+no ', jma2o, jpa2o, jma2o+no |
---|
800 | IF ( jma2o + no .GT. jpa2o ) THEN |
---|
801 | WRITE (nout, *) 'Pas assez de voisins pour la methode essai' |
---|
802 | STOP |
---|
803 | END IF |
---|
804 | !! Liste des points oceans |
---|
805 | no = COUNT ( looce .AND. iomskp == 0_il) |
---|
806 | WRITE (nout,*) 'Nombre de points oce cibles : ', no |
---|
807 | ALLOCATE ( jo_n (no), STAT=ierr) ; CALL chk_allo (ierr, 'jo_n') |
---|
808 | !! |
---|
809 | no = 0 ; z_d_1 = 0.0_rl |
---|
810 | DO jo = 1, jpon |
---|
811 | IF ( .NOT. looce (jo) ) CYCLE |
---|
812 | IF ( iomskp (jo) == 1 ) CYCLE |
---|
813 | no = no + 1 |
---|
814 | jo_n (no) = jo |
---|
815 | z_d_1 = z_d_1 + xosrft (jo) |
---|
816 | END DO |
---|
817 | no = COUNT ( looce .AND. iomskp == 0_il) |
---|
818 | WRITE (nout,*) 'Surface : ', z_d_1 |
---|
819 | SeekAtm_5_1: DO ja = 1, jpan |
---|
820 | IF (itarget (ja) .GE. 1_il) CYCLE SeekAtm_5_1 |
---|
821 | IF (iamskp (ja) .EQ. 1_il) CYCLE SeekAtm_5_1 |
---|
822 | IF (.NOT. laland (ja) ) CYCLE SeekAtm_5_1 |
---|
823 | DO kn = 1, no |
---|
824 | jo = jo_n (kn) |
---|
825 | IF (nvo (jo) < jpa2o ) THEN ! Reste un poids libre |
---|
826 | nvo (jo) = nvo (jo) + 1 |
---|
827 | wa2o (nvo (jo), jo) = xosrft (jo) / z_d_1 |
---|
828 | ka2o (nvo (jo), jo) = ja |
---|
829 | ELSE |
---|
830 | WRITE(nout,*) ' Plus de place ', ja, jot, nvo (jot), z_d_1 |
---|
831 | STOP |
---|
832 | ENDIF |
---|
833 | END DO |
---|
834 | END DO SeekAtm_5_1 |
---|
835 | !! |
---|
836 | CALL gather_wei ('apres essai') |
---|
837 | !! |
---|
838 | CALL verif ( 'Cas essai') |
---|
839 | CALL bilan ("Lessai avant gather ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral") |
---|
840 | CALL bilan ("Lessai avant gather ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral") |
---|
841 | CALL bilan ("Lessai avant gather ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral") |
---|
842 | !! |
---|
843 | WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o |
---|
844 | END IF Essai |
---|
845 | !! |
---|
846 | !! Traitement des points atmosphere en bord de cote |
---|
847 | !! |
---|
848 | !-$$ IfNear: IF (lnear) THEN |
---|
849 | !-$$ WRITE (nout,*) 'Extension de 1 point a l''interieur, vers le point ocean le plus proche (lnear=TRUE)' |
---|
850 | !-$$ SeekAtm1: DO ja = 1, jpan |
---|
851 | !-$$ IF (.NOT. laland_bord(ja) ) CYCLE SeekAtm1 |
---|
852 | !-$$ IF (iamskp(ja) .EQ. 1_il ) CYCLE SeekAtm1 |
---|
853 | !-$$ IF (.NOT. laland(ja) ) CYCLE SeekAtm1 |
---|
854 | !-$$ zzmin = REAL (geodist(0.0_rl, 90.0_rl, 0.0_rl, -90.0_rl), KIND=rl) |
---|
855 | !-$$ ltrouve = .FALSE. ; jot = 0_il |
---|
856 | !-$$ SeekOce1: DO jo = 1, jpon |
---|
857 | !-$$ IF (.NOT. locot(jo) ) CYCLE SeekOce1 |
---|
858 | !-$$ IF (iomskp(jo) .EQ. 1_il ) CYCLE SeekOce1 |
---|
859 | !-$$ z_d_1 = REAL (geodist(xalont(ja), xalatt(ja), xolont(jo), xolatt(jo)), KIND=rl) |
---|
860 | !-$$ IF (z_d_1 .LT. zzmin) THEN |
---|
861 | !-$$ zzmin = z_d_1 |
---|
862 | !-$$ ltrouve = .TRUE. ; jot = jo |
---|
863 | !-$$ END IF |
---|
864 | !-$$ ENDDO SeekOce1 |
---|
865 | !-$$ Trouve11: IF (ltrouve) THEN |
---|
866 | !-$$ WRITE (nout,*) 'Point oce proche ', jot, m2oi(jot), m2oj(jot) |
---|
867 | !-$$ !! Calcul du poids |
---|
868 | !-$$ z_d_1 = xasrft(ja)/xosrft(jot) |
---|
869 | !-$$ !! Rangement des poids |
---|
870 | !-$$ ! On cherche si ce point atm est deja vise par ce point ocean (normalement non) |
---|
871 | !-$$ ltrouve2 = .FALSE. |
---|
872 | !-$$ L11: DO jn = 1, jpa2o |
---|
873 | !-$$ IF (ka2o(jn,jot) == ja ) THEN |
---|
874 | !-$$ wa2o(jn,jot) = wa2o(jn,jot) + z_d_1 |
---|
875 | !-$$ WRITE (nout,*) 'Ajout au point oce ', jo, m2oi(jo), m2oj(jo) |
---|
876 | !-$$ ltrouve2 = .TRUE. |
---|
877 | !-$$ EXIT L11 |
---|
878 | !-$$ ENDIF |
---|
879 | !-$$ ENDDO L11 |
---|
880 | !-$$ Trouve12: IF (.NOT. ltrouve2 ) THEN |
---|
881 | !-$$ IF (nvo(jot) < jpa2o ) THEN ! Reste un poids libre |
---|
882 | !-$$ nvo(jot) = nvo(jot) + 1 |
---|
883 | !-$$ wa2o(nvo(jot),jot) = z_d_1 |
---|
884 | !-$$ ka2o(nvo(jot),jot) = ja |
---|
885 | !-$$ WRITE (nout,*) 'Vers point oce ', jot, m2oi(jot), m2oj(jot) |
---|
886 | !-$$ ELSE ! Trouve2 |
---|
887 | !-$$ WRITE(nout,*) 'Plus de place' |
---|
888 | !-$$ STOP |
---|
889 | !-$$ ENDIF |
---|
890 | !-$$ END IF Trouve12 |
---|
891 | !-$$ ELSE ! Trouve |
---|
892 | !-$$ !! Message d'erreur |
---|
893 | !-$$ WRITE (unit = nout, fmt = '("Pas de voisin proche pour ", 1I5,2I4)' ) ja, m2ai(ja), m2aj(ja) |
---|
894 | !-$$ END IF Trouve11 |
---|
895 | !-$$ END DO SeekAtm1 |
---|
896 | !-$$ !! |
---|
897 | !-$$ CALL verif ( 'Cas Near' ) |
---|
898 | !-$$ !! |
---|
899 | !-$$ WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o |
---|
900 | !-$$ !! |
---|
901 | !-$$ END IF IfNear |
---|
902 | !! |
---|
903 | !-$$ IfNei: IF (lnei) THEN |
---|
904 | !-$$ WRITE(nout,*) 'Route le run-off vers les autres voisins mouilles (lnei=TRUE)' |
---|
905 | !-$$ i2a = 0_il |
---|
906 | !-$$ WHERE (RESHAPE(o2amask,(/jpai,jpaj/)) .GT. 0.0_rl ) i2a = 1_il |
---|
907 | !-$$ SeekAtm2: DO ja = 1_il, jpan |
---|
908 | !-$$ IF (.NOT. laland_bord(ja) ) CYCLE SeekAtm2 |
---|
909 | !-$$ IF (iamskp(ja) .EQ. 1_il ) CYCLE SeekAtm2 |
---|
910 | !-$$ IF (.NOT. laland(ja) ) CYCLE SeekAtm2 |
---|
911 | !-$$ !! Nombre de voisins |
---|
912 | !-$$ jai = m2ai(ja) ; jaj = m2aj(ja) ! Indices 2D |
---|
913 | !-$$ jai_p = mai(jai+1_il) ; jai_m = mai(jai-1_il) |
---|
914 | !-$$ jaj_p = maj(jaj+1_il) ; jaj_m = maj(jaj-1_il) |
---|
915 | !-$$ inum = & |
---|
916 | !-$$ & i2a (jai_p, jaj ) & |
---|
917 | !-$$ & + i2a (jai_m, jaj ) & |
---|
918 | !-$$ & + i2a (jai , jaj_p) & |
---|
919 | !-$$ & + i2a (jai , jaj_m) & |
---|
920 | !-$$ & + i2a (jai_p, jaj_p) & |
---|
921 | !-$$ & + i2a (jai_p, jaj_m) & |
---|
922 | !-$$ & + i2a (jai_m, jaj_p) & |
---|
923 | !-$$ & + i2a (jai_m, jaj_m) |
---|
924 | !-$$ IF (inum == 0) WRITE(nout,*) "erreur ", ja |
---|
925 | !-$$ !! |
---|
926 | !-$$ IF (i2a(jai_p, jaj ) == 1_il ) CALL complet_run (m1a(jai_p, jaj ), m1a(jai,jaj), inum) |
---|
927 | !-$$ IF (i2a(jai_m, jaj ) == 1_il ) CALL complet_run (m1a(jai_m, jaj ), m1a(jai,jaj), inum) |
---|
928 | !-$$ IF (i2a(jai , jaj_p) == 1_il ) CALL complet_run (m1a(jai , jaj_p), m1a(jai,jaj), inum) |
---|
929 | !-$$ IF (i2a(jai , jaj_m) == 1_il ) CALL complet_run (m1a(jai , jaj_m), m1a(jai,jaj), inum) |
---|
930 | !-$$ IF (i2a(jai_p, jaj_p) == 1_il ) CALL complet_run (m1a(jai_p, jaj_p), m1a(jai,jaj), inum) |
---|
931 | !-$$ IF (i2a(jai_p, jaj_m) == 1_il ) CALL complet_run (m1a(jai_p, jaj_m), m1a(jai,jaj), inum) |
---|
932 | !-$$ IF (i2a(jai_m, jaj_p) == 1_il ) CALL complet_run (m1a(jai_m, jaj_p), m1a(jai,jaj), inum) |
---|
933 | !-$$ IF (i2a(jai_m, jaj_m) == 1_il ) CALL complet_run (m1a(jai_m, jaj_m), m1a(jai,jaj), inum) |
---|
934 | !-$$ END DO SeekAtm2 |
---|
935 | !-$$ !! |
---|
936 | !-$$ CALL verif ('Cas Nei') |
---|
937 | !-$$ !! |
---|
938 | !-$$ WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o |
---|
939 | !-$$ !! |
---|
940 | !-$$ END IF IfNei |
---|
941 | !! |
---|
942 | CALL int_wei |
---|
943 | !! |
---|
944 | !! Controle |
---|
945 | !! |
---|
946 | CALL verif (' Normalisation') |
---|
947 | !! |
---|
948 | CALL bilan ("Apres lcoast ", c_case="terre" , c_int_atm="Integral", c_int_oce="Local") |
---|
949 | CALL bilan ("Apres lcoast ", c_case="terre100", c_int_atm="Integral", c_int_oce="Local") |
---|
950 | CALL bilan ("Apres lcoast ", c_case="cote" , c_int_atm="Integral", c_int_oce="Local") |
---|
951 | !! |
---|
952 | WRITE (unit = nout, fmt = * ) 'Nombre de voisins : ', jma2o |
---|
953 | !! |
---|
954 | CALL gather_wei ('apres normalisation') |
---|
955 | CALL verif ('apres gather' ) |
---|
956 | |
---|
957 | !! Atm mask interpolated to oce |
---|
958 | za = REAL (1_il-iamskt, KIND = rl ) |
---|
959 | CALL inter_a2o (za, a2omask ) |
---|
960 | !! 1 on Atm interpolated to oce |
---|
961 | za = REAL (1_il , KIND = rl ) |
---|
962 | CALL inter_a2o (za, a2ofull ) |
---|
963 | |
---|
964 | !! |
---|
965 | !! Checking weights |
---|
966 | WRITE (unit = nout, fmt = '("Neighbors a2o : ", 3I9 )' ) & |
---|
967 | & MINVAL (nvo, nvo /=0), MAXVAL (nvo), jma2o |
---|
968 | WRITE (unit = nout, fmt = '("Index a2o : ", 2I9 )' ) & |
---|
969 | & MINVAL (ka2o, ka2o /= 0), MAXVAL (ka2o) |
---|
970 | jind = MINLOC(wa2o, wa2o /= 0.0_rl) ; jo = jind(2) ; jn = jind(1) |
---|
971 | joi = m2oi(jo) ; joj = m2oj(jo) |
---|
972 | WRITE (unit = nout, fmt = '("Weights a2o mini : ", 4I6, 1E13.4, 3I6)' ) & |
---|
973 | & jo, jn, joi, joj, wa2o(jn,jo), ka2o(jn,jo), m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo)) |
---|
974 | jind = MAXLOC(wa2o, wa2o /= 0.0_rl) ; jo = jind(2) ; jn = jind(1) |
---|
975 | joi = m2oi(jo) ; joj = m2oj(jo) |
---|
976 | WRITE (unit = nout, fmt = '("Weights a2o maxi : ", 4I6, 1E13.4, 3I6)' ) & |
---|
977 | & jo, jn, joi, joj, wa2o(jn,jo), ka2o(jn,jo), m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo)) |
---|
978 | ! |
---|
979 | jind(1:1) = INT(MAXLOC (nvo),il) ; jo = jind(1) |
---|
980 | WRITE (unit = nout, fmt = *) jo, m2oi(jo), m2oj(jo), ka2o(1:nvo(jo),jo), & |
---|
981 | & (m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo)), jn = 1, nvo(jo)) |
---|
982 | !! |
---|
983 | clweight = "WEIGHTS4" ; cladress = "ADRESSE4" |
---|
984 | WRITE (unit = nout, fmt = *) cladress, " ", clweight, " ", jpa2o, jma2o |
---|
985 | !! |
---|
986 | !! Controle |
---|
987 | !! |
---|
988 | WRITE (nout, *) 'Verif avant ecriture' |
---|
989 | WRITE (nout, *) 'wa2o ', MINVAL (wa2o), MAXVAL (wa2o), COUNT (wa2o > 0.0_rl) |
---|
990 | WRITE (nout, *) 'ka2o ', MINVAL (ka2o), MAXVAL (ka2o), COUNT (ka2o > 0 ) |
---|
991 | WRITE (nout, *) 'xasrft ', MINVAL (xasrft), MAXVAL (xasrft) |
---|
992 | WRITE (nout, *) 'o2amask ', MINVAL (o2amask), MAXVAL (o2amask), SUM (o2amask) |
---|
993 | WRITE (nout, *) 'xosrft ', MINVAL (xosrft), MAXVAL (xosrft) |
---|
994 | WRITE (nout, *) 'iomskt ', MINVAL (iomskt), MAXVAL (iomskt), COUNT (iomskt == 0) |
---|
995 | WRITE (nout, *) 'iomskp ', MINVAL (iomskp), MAXVAL (iomskp), COUNT (iomskp == 0) |
---|
996 | WRITE (nout, *) 'iamskt ', MINVAL (iamskt), MAXVAL (iamskt), COUNT (iamskt == 0) |
---|
997 | WRITE (nout, *) 'locot ', COUNT (locot) |
---|
998 | WRITE (nout, *) 'lacot ', COUNT (lacot) |
---|
999 | !! |
---|
1000 | ENDIF L_LCOAST |
---|
1001 | !! |
---|
1002 | END IF E_DRYRUN |
---|
1003 | !! |
---|
1004 | !! |
---|
1005 | WRITE (nout, *) ' ' |
---|
1006 | WRITE (nout, *) ' ' |
---|
1007 | WRITE (nout, *) 'Eventuels traitement specifiques' |
---|
1008 | WRITE (nout, *) ' ' |
---|
1009 | WRITE (nout, *) ' ' |
---|
1010 | |
---|
1011 | |
---|
1012 | !! Extension ?? |
---|
1013 | |
---|
1014 | !! Traitements specifiques, tres specifiques ... !! |
---|
1015 | IF ( TRIM(c_suffix) == "_runoffNord09" ) THEN |
---|
1016 | WRITE (nout, *) 'Traitement specifique _runoffNord09' |
---|
1017 | DO jo = 1, jpon |
---|
1018 | IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.9 |
---|
1019 | END DO |
---|
1020 | END IF |
---|
1021 | ! |
---|
1022 | IF ( TRIM(c_suffix) == "_runoffNord07" ) THEN |
---|
1023 | WRITE (nout, *) 'Traitement specifique _runoffNord07' |
---|
1024 | DO jo = 1, jpon |
---|
1025 | IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.7 |
---|
1026 | END DO |
---|
1027 | END IF |
---|
1028 | |
---|
1029 | IF ( TRIM(c_suffix) == "_runoffNord00" ) THEN |
---|
1030 | WRITE (nout, *) 'Traitement specifique _runoffNord00' |
---|
1031 | DO jo = 1, jpon |
---|
1032 | IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.0 |
---|
1033 | END DO |
---|
1034 | END IF |
---|
1035 | !! |
---|
1036 | !! |
---|
1037 | !! Verif nombre de voisins |
---|
1038 | !! |
---|
1039 | IF ( jma2or /= jma2o ) THEN |
---|
1040 | WRITE (nout, *) 'Erreur nombre de voisins dans run.def : ' |
---|
1041 | WRITE (nout, *) ' jma2o calcule / lu : ', jma2o, jma2or |
---|
1042 | !STOP |
---|
1043 | ENDIF |
---|
1044 | !! |
---|
1045 | !! Ecriture des poids |
---|
1046 | !! |
---|
1047 | WRITE (unit = c_r, fmt = '(1i1)' ) rk_out |
---|
1048 | ! |
---|
1049 | IF ( l_wei_i4) THEN |
---|
1050 | WRITE (unit = c_i, fmt = '(1i1)' ) i_4 |
---|
1051 | cfname4 = TRIM(mo_name) // ".wa2o.runoff.i" // TRIM(c_i) // ".r" // TRIM(c_r) // TRIM(c_suffix) |
---|
1052 | OPEN (unit = nwei4a2o, file = TRIM(cfname4), form = 'unformatted', status = 'replace',& |
---|
1053 | & action = 'write') |
---|
1054 | WRITE (nout, *) 'Poids i_4 a -> o: ', TRIM(cfname4) |
---|
1055 | END IF |
---|
1056 | ! |
---|
1057 | IF (l_wei_i8) THEN |
---|
1058 | WRITE (unit = c_i, fmt = '(1i1)' ) i_8 |
---|
1059 | cfname8 = TRIM(mo_name) // ".wa2o.runoff.i" // TRIM(c_i) // ".r" // TRIM(c_r) // TRIM(c_suffix) |
---|
1060 | OPEN (unit = nwei8a2o, file = TRIM(cfname8), form = 'unformatted', status = 'replace',& |
---|
1061 | & action = 'write') |
---|
1062 | WRITE (nout, *) 'Poids i_8 a -> o : ', TRIM(cfname8) |
---|
1063 | END IF |
---|
1064 | ! |
---|
1065 | cldiag_a2o = TRIM(a2o_name) // '.runcoa.diag' |
---|
1066 | clw_a2o = TRIM(mo_name) // '.wa2o.runcoa' |
---|
1067 | clw_a2o_mct= TRIM(mo_name) // '.wa2o_mct.runcoa' |
---|
1068 | ! |
---|
1069 | CALL wri_weights_a2o (cldiag_a2o, clw_a2o, clw_a2o_mct, "4", l_fulldiag, lo_src_grid_frac=.FALSE., lo_dst_grid_frac=.FALSE., & |
---|
1070 | & co_omsk=TRIM(cotes_omsk), co_amsk=TRIM(cotes_amsk) ) |
---|
1071 | ! |
---|
1072 | IF (lriv) THEN |
---|
1073 | !! Dataset 3 : cas particulier du run off des rivieres |
---|
1074 | !! |
---|
1075 | ka2o = 0 ; wa2o = 0.0e0_rl ; nvo = 0 ; nva = 0 ; jma2o = 1 |
---|
1076 | OPEN (unit = nriv, file = '/home/clim/ocean/CPL/DATA/riviere.dat', status = 'old', action = 'read', & |
---|
1077 | & form = 'formatted', position = 'rewind' ) |
---|
1078 | jovoid = 1 |
---|
1079 | DO jr = 1, jpr |
---|
1080 | READ (unit = nriv, fmt = '(1A1,2I3,1I5,1A30)' ) clriv, joi, joj, nlmd, clname |
---|
1081 | clname = ADJUSTL (TRIM (ADJUSTL (clname))) |
---|
1082 | jai = m2ai (nlmd) ; jaj = m2aj (nlmd) |
---|
1083 | jaj = jpaj - jaj + 1 |
---|
1084 | ja = m1a(jai,jaj) |
---|
1085 | IF (joi == 99 .AND. joj == 99 ) THEN |
---|
1086 | !! On ramene sur la ligne j=1 de l'ocean pour faire un bilan plus tard |
---|
1087 | jovoid = jovoid + 1 |
---|
1088 | WRITE (unit = nout, fmt = & |
---|
1089 | & '(1A, " (", 1A30, ") ", 1I3, " : ", 4I5, " : ", 2F6.1 )' ) & |
---|
1090 | & clriv, clname, jr, nlmd, ja, jai, jaj, xalont (ja), xalatt (ja) |
---|
1091 | joi = jovoid ; joj = 1 |
---|
1092 | jo = m1o (joi, joj) |
---|
1093 | nva (ja) = nva (ja) + 1 |
---|
1094 | nvo (jo) = nvo (jo) + 1 |
---|
1095 | ka2o (nvo (jo), jo) = ja |
---|
1096 | wa2o (nvo (jo), jo) = 1.0e0_rl |
---|
1097 | ELSE |
---|
1098 | jo = m1o (joi, joj) |
---|
1099 | nva (ja) = nva (ja) + 1 |
---|
1100 | nvo (jo) = nvo (jo) + 1 |
---|
1101 | ka2o (nvo (jo), jo) = ja |
---|
1102 | IF (lint_atm .AND. lint_oce ) THEN |
---|
1103 | !! Atm envoie : kg.s-1, Oce reçoit : kg.s-1 |
---|
1104 | wa2o(nvo(jo),jo) = 1.0e0_rl |
---|
1105 | ELSE |
---|
1106 | !! Atm envoie kg/m ^2/s. Oce reçoit kg/m^2/s |
---|
1107 | wa2o (nvo (jo), jo) = xasrft (ja) / xosrft (jo) |
---|
1108 | END IF |
---|
1109 | WRITE (unit = nout, fmt = & |
---|
1110 | & '(1A, " (", 1A30, ") ", 1I3, " : ", 4I5, " -> ", 2I4, 1I2, " : ", 1E12.4, " ", 2F6.1, " -> ", 2F6.1)' ) & |
---|
1111 | & clriv, clname, jr, nlmd, ja, jai, jaj, joi, joj, nvo(jo), wa2o (nvo(jo),jo), xalont(ja), xalatt(ja), & |
---|
1112 | & clo_lon(xolont(jo),xalont(ja)), xolatt (jo) |
---|
1113 | END IF |
---|
1114 | END DO |
---|
1115 | CLOSE (unit = nriv) |
---|
1116 | jma2o = MAXVAL (nvo) |
---|
1117 | WRITE (unit = nout, fmt = '("Voisins atm : ", 2I9 )' ) & |
---|
1118 | & MINVAL(nva), MAXVAL(nva) |
---|
1119 | WRITE (unit = nout, fmt = '("Voisins ocean : ", 2I9 )' ) & |
---|
1120 | & MINVAL(nvo), MAXVAL(nvo) |
---|
1121 | WRITE (unit = nout, fmt = '("Point bidons : ", 2I9 )' ) jovoid |
---|
1122 | !! |
---|
1123 | END IF |
---|
1124 | |
---|
1125 | !! |
---|
1126 | !! Verif nombre de voisins |
---|
1127 | !! |
---|
1128 | IF ( jma2or /= jma2o ) THEN |
---|
1129 | WRITE (nout, *) 'Erreur nombre de voisins dans run.def : ' |
---|
1130 | WRITE (nout, *) ' jma2o calcule / lu : ', jma2o, jma2or |
---|
1131 | !STOP |
---|
1132 | ENDIF |
---|
1133 | !! |
---|
1134 | !! Ecriture des poids |
---|
1135 | !! |
---|
1136 | !! Ouverture des fichiers poids |
---|
1137 | !! |
---|
1138 | WRITE (unit = c_r, fmt = '(1i1)' ) rk_out |
---|
1139 | ! |
---|
1140 | IF (l_wei_i4) THEN |
---|
1141 | WRITE (unit = c_i, fmt = '(1i1)' ) i_4 |
---|
1142 | cfname4 = TRIM(mo_name) // ".wa2o.rivflu" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r) |
---|
1143 | OPEN (unit = nwei4a2o, file = TRIM(cfname4), form = 'unformatted', status = 'replace',& |
---|
1144 | & action = 'write') |
---|
1145 | WRITE (nout, *) 'Poids i4 a -> o: ', TRIM(cfname4) |
---|
1146 | END IF |
---|
1147 | ! |
---|
1148 | IF (l_wei_i4) THEN |
---|
1149 | WRITE (unit = c_i, fmt = '(1i1)' ) i_8 |
---|
1150 | cfname8 = TRIM (mo_name) // ".wa2o.rivflu" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r) |
---|
1151 | OPEN (unit = nwei8a2o, file = TRIM(cfname8), form = 'unformatted', status = 'replace',& |
---|
1152 | & action = 'write') |
---|
1153 | WRITE (nout, *) 'Poids i8 a -> o : ', TRIM(cfname8) |
---|
1154 | END IF |
---|
1155 | |
---|
1156 | ! |
---|
1157 | cldiag_a2o = TRIM (a2o_name) // '.rivflu.diag' |
---|
1158 | clw_a2o = TRIM (mo_name) // '.wa2o.rivflu' |
---|
1159 | clw_a2o_mct = "rmp_"//TRIM(camod_t)//"_to_"//TRIM(comod_t)//"_MOSAIC_rivflu" |
---|
1160 | ! |
---|
1161 | CALL wri_weights_a2o (cldiag_a2o, clw_a2o, clw_a2o_mct, "3", l_fulldiag, lo_src_grid_frac=.FALSE., lo_dst_grid_frac=.FALSE., & |
---|
1162 | & co_omsk='noperio', co_amsk='full' ) |
---|
1163 | |
---|
1164 | !! |
---|
1165 | IF (l_wei_i4) CLOSE (unit = nwei8a2o) |
---|
1166 | IF (l_wei_i8) CLOSE (unit = nwei4a2o) |
---|
1167 | !! |
---|
1168 | WRITE (nout,*) 'Appel flioclo' |
---|
1169 | CALL flioclo |
---|
1170 | !! |
---|
1171 | STOP |
---|
1172 | !! |
---|
1173 | CONTAINS |
---|
1174 | !! |
---|
1175 | SUBROUTINE diag_itarget (clmess, lwrite) |
---|
1176 | CHARACTER (len=*), INTENT (in) :: clmess |
---|
1177 | LOGICAL, INTENT (in) :: lwrite |
---|
1178 | INTEGER, SAVE :: nappel = 0 |
---|
1179 | CHARACTER (len=2) :: cc |
---|
1180 | !! |
---|
1181 | nappel = nappel + 1 |
---|
1182 | !! |
---|
1183 | itarget = 0_il ; wtarget = 0.0_rl |
---|
1184 | DO jo = 1, jpon |
---|
1185 | DO jn = 1, jpa2o |
---|
1186 | ja = ka2o (jn, jo) |
---|
1187 | IF (ja /= 0 ) THEN |
---|
1188 | itarget (ja) = itarget (ja) + 1_il |
---|
1189 | wtarget (ja) = wtarget (ja) + wa2o (jn, jo) |
---|
1190 | ENDIF |
---|
1191 | END DO |
---|
1192 | END DO |
---|
1193 | !! |
---|
1194 | !CALL lbc (itarget, naperio, 'T', jpai, jpaj) |
---|
1195 | !CALL lbc (wtarget, naperio, 'T', jpai, jpaj) |
---|
1196 | !! |
---|
1197 | IF (lwrite) THEN |
---|
1198 | !$$$ WRITE(nout,*) 'itarget, etape '// clmess |
---|
1199 | !$$$ CALL prihin (RESHAPE(itarget,(/jpai,jpaj/)),1_il,nout) |
---|
1200 | WRITE (unit = cc, fmt = '(1I2.2)' ) nappel |
---|
1201 | CALL fliodefv (nid_it, 'itarget' // cc, (/1, 2/), v_t=flio_i4) |
---|
1202 | CALL fliodefv (nid_it, 'ntarget' // cc, (/1, 2/), v_t=flio_i4) |
---|
1203 | CALL fliodefv (nid_it, 'wtarget' // cc, (/1, 2/), v_t=flio_r4) |
---|
1204 | CALL fliodefv (nid_it, 'non_target' // cc, (/1, 2/), v_t=flio_i4) |
---|
1205 | |
---|
1206 | |
---|
1207 | CALL flioputa (nid_it, 'itarget' // cc, 'long_name ', clmess // ' itarget' // cc) |
---|
1208 | CALL flioputv (nid_it, 'itarget' // cc, MIN (1, RESHAPE (itarget, (/jpai, jpaj/))) ) |
---|
1209 | |
---|
1210 | CALL flioputa (nid_it, 'ntarget' // cc, 'long_name ', clmess // ' ntarget' // cc) |
---|
1211 | CALL flioputv (nid_it, 'ntarget' // cc, RESHAPE (itarget, (/jpai, jpaj/)) ) |
---|
1212 | |
---|
1213 | CALL flioputa (nid_it, 'wtarget' // cc, 'long_name ', clmess // ' wtarget' // cc) |
---|
1214 | CALL flioputv (nid_it, 'wtarget' // cc, RESHAPE (wtarget, (/jpai, jpaj/)) ) |
---|
1215 | |
---|
1216 | wtarget = 0 |
---|
1217 | WHERE (itarget .EQ. 0 .AND. (laland .OR. lacot)) wtarget = 1.0_rl |
---|
1218 | CALL flioputa (nid_it, 'non_target' // cc, 'long_name ', clmess // ' non_target' // cc) |
---|
1219 | CALL flioputv (nid_it, 'non_target' //cc, RESHAPE (NINT (wtarget), (/jpai, jpaj/)) ) |
---|
1220 | END IF |
---|
1221 | !! |
---|
1222 | END SUBROUTINE diag_itarget |
---|
1223 | !! |
---|
1224 | SUBROUTINE int_wei |
---|
1225 | !! |
---|
1226 | IMPLICIT NONE |
---|
1227 | WRITE(nout,*) 'Debut normalisation' |
---|
1228 | !! |
---|
1229 | IF (.NOT. lint_atm) THEN |
---|
1230 | WRITE(nout,*) 'Multiplication par les surfaces atm' |
---|
1231 | DO jo = 1, jpon |
---|
1232 | DO jn = 1, jma2o |
---|
1233 | ja = ka2o (jn, jo) |
---|
1234 | IF (ja > 0 ) wa2o (jn,jo) = wa2o (jn,jo) * xasrft (ja) |
---|
1235 | END DO |
---|
1236 | ENDDO |
---|
1237 | ENDIF |
---|
1238 | !! |
---|
1239 | IF (.NOT. lint_oce) THEN |
---|
1240 | WRITE (nout,*) 'Division par les surfaces oce' |
---|
1241 | DO jo = 1, jpon |
---|
1242 | DO jn = 1, jma2o |
---|
1243 | wa2o (jn,jo) = wa2o (jn,jo) / xosrft (jo) |
---|
1244 | END DO |
---|
1245 | ENDDO |
---|
1246 | ENDIF |
---|
1247 | END SUBROUTINE int_wei |
---|
1248 | !! |
---|
1249 | SUBROUTINE complet_run (ka1, ka2, kn) |
---|
1250 | INTEGER (kind=il), INTENT(in) :: ka1, ka2 ! Index of atm boxes |
---|
1251 | INTEGER (kind=il), INTENT(in) :: kn ! Number of neighbors |
---|
1252 | INTEGER (kind=il) :: jt |
---|
1253 | !! |
---|
1254 | DO jo = 1, jpon |
---|
1255 | IF (iomskt(jo) == 1 ) CYCLE |
---|
1256 | IF (iomskp(jo) .EQ. 1_il ) CYCLE |
---|
1257 | DO jn = 1, jpa2o |
---|
1258 | IF (ka2o(jn,jo) == ka1 ) THEN |
---|
1259 | nvo(jo) = nvo(jo)+1 |
---|
1260 | IF (nvo(jo) > jpa2o ) THEN |
---|
1261 | WRITE(nout,*) 'Pas assez de voisins pour la completion ' |
---|
1262 | WRITE(nout,fmt='(1I7,2I4,2F10.1)' ) ka1, m2ai(ka1), m2aj(ka1), & |
---|
1263 | & xalont(ka1), xalatt(ka1) |
---|
1264 | WRITE(nout,fmt='(1I7,2I4,2F10.1)' ) ka2, m2ai(ka2), m2aj(ka2), & |
---|
1265 | & xalont(ka2), xalatt(ka2) |
---|
1266 | WRITE(nout,fmt='(1I7,2I4,2F10.1,2I8)') jo, m2oi(jo), m2oj(jo), & |
---|
1267 | & xolont(jo), xolatt(jo), jn, nvo(jo) |
---|
1268 | DO jt = 1, jpa2o |
---|
1269 | WRITE(nout,fmt=*) jt, ka2o(jt,jo), xalont(ka2o(jt,jo)), xalatt(ka2o(jt,jo)), & |
---|
1270 | & lacot(ka2o(jt,jo)) |
---|
1271 | ENDDO |
---|
1272 | END IF |
---|
1273 | ka2o (nvo (jo), jo) = ka2 |
---|
1274 | wa2o (nvo (jo), jo) = wa2o(jn,jo) * xasrft(ka2) / xasrft(ka1) / REAL(kn,KIND=rl) |
---|
1275 | END IF |
---|
1276 | END DO |
---|
1277 | ENDDO |
---|
1278 | END SUBROUTINE complet_run |
---|
1279 | !! |
---|
1280 | SUBROUTINE gather_wei (clmess) |
---|
1281 | CHARACTER (len=*), INTENT (in) :: clmess |
---|
1282 | INTEGER (kind=il) :: jn1, jn2, jo2 |
---|
1283 | !! Eliminate redundancy |
---|
1284 | DO jo = 1, jpon |
---|
1285 | DO jn1 = 1, jpa2o-1 |
---|
1286 | IF (ka2o (jn1,jo) == 0 ) CYCLE |
---|
1287 | DO jn2 = jn1+1, jpa2o |
---|
1288 | IF (ka2o (jn1, jo) == ka2o (jn2, jo) ) THEN |
---|
1289 | WRITE(nout,*) 'Redondance ', jo, jn1, jn2, ka2o(jn1,jo), ka2o (jn2,jo), wa2o(jn1,jo), wa2o(jn2,jo) |
---|
1290 | wa2o (jn1, jo) = wa2o (jn1, jo) + wa2o (jn2, jo) |
---|
1291 | ka2o (jn2, jo) = 0 |
---|
1292 | END IF |
---|
1293 | END DO |
---|
1294 | END DO |
---|
1295 | END DO |
---|
1296 | !! |
---|
1297 | !! Seek for strangeness |
---|
1298 | DO jo = 1, jpon |
---|
1299 | DO jn = 1, jpa2o |
---|
1300 | ja = ka2o (jn, jo) |
---|
1301 | IF (ja == 0) CYCLE |
---|
1302 | IF (wa2o (jn, jo) == 0.0E0 ) THEN |
---|
1303 | WRITE (nout, *) 'Null weight : ', jo, jn, ja, locot (jo), lacot (ja) |
---|
1304 | END IF |
---|
1305 | END DO |
---|
1306 | END DO |
---|
1307 | !! |
---|
1308 | !! Periodicity |
---|
1309 | !DO jn = 1, jpa2o |
---|
1310 | ! CALL lbc (wa2o (jn,:), noperio, 'T', jpoi, jpoj) |
---|
1311 | ! CALL lbc (ka2o (jn,:), noperio, 'T', jpoi, jpoj) |
---|
1312 | !END DO |
---|
1313 | CALL diag_itarget (' gather/perio - ' // clmess, lwrite=.TRUE.) |
---|
1314 | !! |
---|
1315 | !! Gather weights |
---|
1316 | DO jo = 1, jpon |
---|
1317 | ktemp = ka2o (:,jo) ; wtemp = wa2o (:,jo) |
---|
1318 | ka2o (:,jo) = 0_il ; wa2o (:,jo) = 0.0_rl |
---|
1319 | kn = 1_il |
---|
1320 | DO jo2 = 1, jpa2o |
---|
1321 | IF (ktemp (jo2) /= 0_il ) THEN |
---|
1322 | ka2o (kn, jo) = ktemp (jo2) |
---|
1323 | wa2o (kn, jo) = wtemp (jo2) |
---|
1324 | kn = kn + 1_il |
---|
1325 | END IF |
---|
1326 | END DO |
---|
1327 | END DO |
---|
1328 | nvo = COUNT (ka2o /= 0_il, DIM = 1) |
---|
1329 | jma2o = MAXVAL (nvo) |
---|
1330 | ! |
---|
1331 | !! |
---|
1332 | !! Periodicity |
---|
1333 | !DO jn = 1, jpa2o |
---|
1334 | ! CALL lbc (wa2o(jn,:), noperio, 'T', jpoi, jpoj) |
---|
1335 | ! CALL lbc (ka2o(jn,:), noperio, 'T', jpoi, jpoj) |
---|
1336 | !END DO |
---|
1337 | !! |
---|
1338 | nvo = COUNT (ka2o /= 0_il, DIM = 1) |
---|
1339 | jma2o = MAXVAL (nvo) |
---|
1340 | !! |
---|
1341 | END SUBROUTINE gather_wei |
---|
1342 | !! |
---|
1343 | SUBROUTINE verif (clmess, l_non_a) |
---|
1344 | !! |
---|
1345 | CHARACTER (len=*), INTENT (in) :: clmess |
---|
1346 | INTEGER :: ja |
---|
1347 | LOGICAL, INTENT (in), OPTIONAL :: l_non_a |
---|
1348 | !! |
---|
1349 | nvo = COUNT (ka2o /= 0_il, DIM = 1) |
---|
1350 | jma2o = MAXVAL (nvo) |
---|
1351 | CALL diag_itarget (clmess, lwrite=.TRUE.) |
---|
1352 | !! |
---|
1353 | WRITE (nout, *) ' --- Verifs -- : ', TRIM (clmess) |
---|
1354 | WRITE (nout, *) 'wa2o ', MINVAL (wa2o), MAXVAL (wa2o), COUNT (wa2o > 0.0_rl) |
---|
1355 | WRITE (nout, *) 'ka2o ', MINVAL (ka2o), MAXVAL (ka2o), COUNT (ka2o > 0 ) |
---|
1356 | WRITE (nout, *) 'xasrft ', MINVAL (xasrft), MAXVAL (xasrft) |
---|
1357 | WRITE (nout, *) 'o2amask ', MINVAL (o2amask), MAXVAL (o2amask), SUM (o2amask) |
---|
1358 | WRITE (nout, *) 'xosrft ', MINVAL (xosrft), MAXVAL (xosrft) |
---|
1359 | WRITE (nout, *) 'iomskt ', MINVAL (iomskt), MAXVAL (iomskt), COUNT (iomskt == 0) |
---|
1360 | WRITE (nout, *) 'iomskp ', MINVAL (iomskp), MAXVAL (iomskp), COUNT (iomskp == 0) |
---|
1361 | WRITE (nout, *) 'iamskt ', MINVAL (iamskt), MAXVAL (iamskt), COUNT (iamskt == 0) |
---|
1362 | WRITE (nout, *) 'locot ', COUNT (locot) |
---|
1363 | WRITE (nout, *) 'lacot ', COUNT (lacot) |
---|
1364 | WRITE (nout, *) 'Points non attribues : ', COUNT (itarget .EQ. 0 .AND. (laland .OR. lacot)) |
---|
1365 | |
---|
1366 | IF (PRESENT (l_non_a)) THEN |
---|
1367 | DO ja = 1, jpan |
---|
1368 | IF (itarget (ja) == 0 .AND. iamskp(ja) == 0 .AND. (laland (ja) .OR. lacot(ja) )) THEN |
---|
1369 | !-$$ WRITE (UNIT=nout,FMT='("Point atm non attribues : ", 3I7, 2F7.1, 1I4)' ) ja, m2ai(ja), m2aj (ja), xalont(ja), xalatt(ja), itarget(ja) |
---|
1370 | WRITE (nout,*) ja, m2ai(ja), m2aj (ja), xalont(ja), xalatt(ja), itarget(ja), laland(ja), lacot(ja), laoce(ja) |
---|
1371 | ENDIF |
---|
1372 | END DO |
---|
1373 | END IF |
---|
1374 | END SUBROUTINE verif |
---|
1375 | !! |
---|
1376 | !! |
---|
1377 | SUBROUTINE noroute |
---|
1378 | !! Determine les points atmospheres non routes |
---|
1379 | lanoroute = .TRUE. |
---|
1380 | |
---|
1381 | DO jn = 1, jpa2o |
---|
1382 | DO jo = 1, jpon |
---|
1383 | ja = ka2o ( jo, jn) |
---|
1384 | IF ( ja /= 0 .AND. wa2o (jo,jn) > 0.0_rl ) THEN |
---|
1385 | lanoroute = .FALSE. |
---|
1386 | END IF |
---|
1387 | END DO |
---|
1388 | END DO |
---|
1389 | |
---|
1390 | END SUBROUTINE noroute |
---|
1391 | !! |
---|
1392 | SUBROUTINE bilan (clmess, c_case, c_int_atm, c_int_oce) |
---|
1393 | !! |
---|
1394 | CHARACTER (len=*), INTENT (in) :: clmess |
---|
1395 | CHARACTER (len=*), INTENT (in) :: c_case ! Type de champ d'entree |
---|
1396 | CHARACTER (len=*), INTENT (in) :: c_int_atm ! Type d'integrale sur l'atm |
---|
1397 | CHARACTER (len=*), INTENT (in) :: c_int_oce ! Type d'integrale sur l'ocean |
---|
1398 | !! |
---|
1399 | WRITE (nout,*) "Bilan :" // TRIM (clmess) // ", cas :", TRIM (c_case), ", c_int_atm :", TRIM (c_int_atm),& |
---|
1400 | & ", c_int_oce :", TRIM(c_int_oce), "." |
---|
1401 | !! |
---|
1402 | za = 0.0_rl |
---|
1403 | !! |
---|
1404 | SELECT CASE (TRIM (c_case)) |
---|
1405 | CASE ('complet') |
---|
1406 | za = 1.0_rl * REAL (1-iamskp, KIND=rl) |
---|
1407 | SELECT CASE (TRIM (c_int_atm)) |
---|
1408 | CASE ('Integral') |
---|
1409 | za = za * xasrft |
---|
1410 | zsuma = SUM (za) |
---|
1411 | CASE ('Local') |
---|
1412 | zsuma = DOT_PRODUCT (za, xasrft) |
---|
1413 | CASE default |
---|
1414 | WRITE (nout,*) 'Erreur type de bilan atm' ; STOP |
---|
1415 | END SELECT |
---|
1416 | !! |
---|
1417 | CASE ('terre') |
---|
1418 | WHERE (laland .OR. lacot) |
---|
1419 | za = 1.0_rl * REAL (1-iamskp, KIND=rl) |
---|
1420 | END WHERE |
---|
1421 | SELECT CASE (TRIM (c_int_atm)) |
---|
1422 | CASE ('Integral') |
---|
1423 | za = za * xasrft |
---|
1424 | zsuma = SUM (za) |
---|
1425 | CASE ('Local') |
---|
1426 | zsuma = DOT_PRODUCT (za, xasrft*(1.0-o2amask)) |
---|
1427 | CASE default |
---|
1428 | WRITE (nout,*) 'Erreur type de bilan atm' ; STOP |
---|
1429 | END SELECT |
---|
1430 | !! |
---|
1431 | CASE ('terre100') |
---|
1432 | WHERE (laland) |
---|
1433 | za = 1.0_rl * REAL (1-iamskp, KIND=rl) |
---|
1434 | END WHERE |
---|
1435 | SELECT CASE (TRIM (c_int_atm)) |
---|
1436 | CASE ('Integral') |
---|
1437 | za = za * xasrft |
---|
1438 | zsuma = SUM (za) |
---|
1439 | CASE ('Local') |
---|
1440 | zsuma = DOT_PRODUCT (za, xasrft) |
---|
1441 | CASE default |
---|
1442 | WRITE (nout,*) 'Erreur type de bilan atm' ; STOP |
---|
1443 | END SELECT |
---|
1444 | !! |
---|
1445 | CASE ('ocean') |
---|
1446 | WHERE (lacot) |
---|
1447 | za = 1.0_rl * REAL (1-iamskp, KIND=rl) |
---|
1448 | END WHERE |
---|
1449 | SELECT CASE (TRIM (c_int_atm)) |
---|
1450 | CASE ('Integral') |
---|
1451 | za = za * xasrft |
---|
1452 | zsuma = SUM (za) |
---|
1453 | CASE ('Local') |
---|
1454 | zsuma = DOT_PRODUCT (za, xasrft*o2amask) |
---|
1455 | CASE default |
---|
1456 | WRITE (nout,*) TRIM(clmess), ' : Erreur type de bilan atm' ; STOP |
---|
1457 | END SELECT |
---|
1458 | !! |
---|
1459 | CASE ('cote') |
---|
1460 | WHERE (lacot) |
---|
1461 | za = 1.0_rl * REAL (1-iamskp, KIND=rl) |
---|
1462 | END WHERE |
---|
1463 | SELECT CASE (TRIM (c_int_atm)) |
---|
1464 | CASE ('Integral') |
---|
1465 | za = za * xasrft |
---|
1466 | zsuma = SUM (za) |
---|
1467 | CASE ('Local') |
---|
1468 | zsuma = DOT_PRODUCT (za, xasrft*o2amask) |
---|
1469 | CASE default |
---|
1470 | WRITE (nout,*) 'Erreur type de bilan atm' ; STOP |
---|
1471 | END SELECT |
---|
1472 | !! |
---|
1473 | CASE default |
---|
1474 | WRITE (nout,*)'Erreur choix de cas' ; STOP |
---|
1475 | STOP |
---|
1476 | END SELECT |
---|
1477 | !! |
---|
1478 | CALL inter_a2o (za, zo) |
---|
1479 | !! |
---|
1480 | SELECT CASE (TRIM (c_int_oce)) |
---|
1481 | CASE ('Integral') |
---|
1482 | zsumo = DOT_PRODUCT (zo, REAL(1-iomskp, KIND=rl)) |
---|
1483 | CASE ('Local') |
---|
1484 | zsumo = DOT_PRODUCT (zo*xosrft,REAL(1-iomskp, KIND=rl)) |
---|
1485 | CASE default |
---|
1486 | WRITE (nout,*) TRIM (clmess), ' : Erreur type de bilan ' ; STOP |
---|
1487 | END SELECT |
---|
1488 | !! |
---|
1489 | WRITE (nout, fmt='(" bilan atm: ", 1E15.6, ", bilan oce: ", 1E15.6)' ) zsuma, zsumo |
---|
1490 | !! |
---|
1491 | RETURN |
---|
1492 | !! |
---|
1493 | END SUBROUTINE bilan |
---|
1494 | !! |
---|
1495 | END PROGRAM cote |
---|
1496 | |
---|