[3326] | 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 | |
---|