Changeset 1319 for LMDZ4/trunk/libf/dyn3dpar/startvar.F90
- Timestamp:
- Feb 23, 2010, 10:29:54 PM (14 years ago)
- File:
-
- 1 moved
-
LMDZ4/trunk/libf/dyn3dpar/startvar.F90 (moved) (moved from LMDZ4/trunk/libf/dyn3dpar/startvar.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/trunk/libf/dyn3dpar/startvar.F90
r1318 r1319 2 2 ! $Id$ 3 3 ! 4 MODULE startvar 4 !******************************************************************************* 5 ! 6 MODULE startvar 7 ! 8 !******************************************************************************* 9 ! 10 !------------------------------------------------------------------------------- 11 ! Purpose: Access data from the database of atmospheric to initialize the model. 12 !------------------------------------------------------------------------------- 13 ! Comments: 14 ! 15 ! * This module is designed to work for Earth (and with ioipsl) 16 ! 17 ! * There are three ways to acces data, depending on the type of field 18 ! which needs to be extracted. In any case the call should come after a restget 19 ! and should be of the type : CALL startget(...) 20 ! 21 ! - A 1D variable on the physical grid : 22 ! CALL startget(varname, iml, jml, lon_in, lat_in, nbindex, & 23 ! champ, val_exp, jml2, lon_in2, lat_in2, ibar ) 24 ! 25 ! - A 2D variable on the dynamical grid : 26 ! CALL startget(varname, iml, jml, lon_in, lat_in, & 27 ! champ, val_exp, jml2, lon_in2, lat_in2, ibar ) 28 ! 29 ! - A 3D variable on the dynamical grid : 30 ! CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, & 31 ! champ, val_exp, jml2, lon_in2, lat_in2, ibar ) 32 ! 33 ! * Data needs to be in NetCDF format 34 ! 35 ! * Variables should have the following names in the files: 36 ! 'RELIEF' : High resolution orography 37 ! 'ST' : Surface temperature 38 ! 'CDSW' : Soil moisture 39 ! 'Z' : Surface geopotential 40 ! 'SP' : Surface pressure 41 ! 'U' : East ward wind 42 ! 'V' : Northward wind 43 ! 'TEMP' : Temperature 44 ! 'R' : Relative humidity 45 ! 46 ! * There is a big mess with the longitude size. Should it be iml or iml+1 ? 47 ! I have chosen to use the iml+1 as an argument to this routine and we declare 48 ! internaly smaller fields when needed. This needs to be cleared once and for 49 ! all in LMDZ. A convention is required. 50 !------------------------------------------------------------------------------- 5 51 #ifdef CPP_EARTH 6 ! This module is designed to work for Earth (and with ioipsl) 7 ! 8 ! 9 ! There are three ways to access data from the database of atmospheric data which 10 ! can be used to initialize the model. This depends on the type of field which needs 11 ! to be extracted. In any case the call should come after a restget and should be of the type : 12 ! CALL startget(...) 13 ! 14 ! We will details the possible arguments to startget here : 15 ! 16 ! - A 2D variable on the dynamical grid : 17 ! CALL startget(varname, iml, jml, lon_in, lat_in, champ, val_ex, jml2, lon_in2, lat_in2, interbar ) 18 ! 19 ! - A 1D variable on the physical grid : 20 ! CALL startget(varname, iml, jml, lon_in, lat_in, nbindex, champ, val_exp, jml2, lon_in2, lat_in2, interbar ) 21 ! 22 ! 23 ! - A 3D variable on the dynamical grid : 24 ! CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, champ, val_exp, jml2, lon_in2, lat_in2, interbar ) 25 ! 26 ! 27 ! There is special constraint on the atmospheric data base except that the 28 ! the data needs to be in netCDF and the variables should have the the following 29 ! names in the file : 30 ! 31 ! 'RELIEF' : High resolution orography 32 ! 'ST' : Surface temperature 33 ! 'CDSW' : Soil moisture 34 ! 'Z' : Surface geopotential 35 ! 'SP' : Surface pressure 36 ! 'U' : East ward wind 37 ! 'V' : Northward wind 38 ! 'TEMP' : Temperature 39 ! 'R' : Relative humidity 40 ! 41 USE ioipsl 42 ! 43 ! 44 IMPLICIT NONE 45 ! 46 ! 47 PRIVATE 48 PUBLIC startget 49 ! 50 ! 51 INTERFACE startget 52 MODULE PROCEDURE startget_phys2d, startget_phys1d, startget_dyn 53 END INTERFACE 54 ! 55 INTEGER, SAVE :: fid_phys, fid_dyn 56 INTEGER, SAVE :: iml_phys, iml_rel, iml_dyn 57 INTEGER, SAVE :: jml_phys, jml_rel, jml_dyn 58 INTEGER, SAVE :: llm_dyn, ttm_dyn 59 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: lon_phys, lon_rug, 60 . lon_alb, lon_rel, lon_dyn 61 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: lat_phys, lat_rug, 62 . lat_alb, lat_rel, lat_dyn 63 REAL, ALLOCATABLE, SAVE, DIMENSION (:) :: levdyn_ini 64 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: relief, zstd, zsig, 65 . zgam, zthe, zpic, zval 66 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rugo, masque, phis 67 ! 68 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: tsol, qsol, psol_dyn 69 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: var_ana3d 70 ! 71 CONTAINS 72 ! 73 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 74 ! 75 SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in, 76 . champ, val_exp, jml2, lon_in2, lat_in2 , interbar, masque_lu ) 77 ! 78 ! There is a big mess with the size in logitude, should it be iml or iml+1. 79 ! I have chosen to use the iml+1 as an argument to this routine and we declare 80 ! internaly smaler fields when needed. This needs to be cleared once and for all in LMDZ. 81 ! A convention is required. 82 ! 83 ! 84 CHARACTER*(*), INTENT(in) :: varname 85 INTEGER, INTENT(in) :: iml, jml ,jml2 86 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 87 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 88 REAL, INTENT(inout) :: champ(iml,jml) 89 REAL, INTENT(in) :: val_exp 90 REAL, INTENT(in), optional :: masque_lu(iml,jml) 91 LOGICAL interbar 92 ! 93 ! This routine only works if the variable does not exist or is constant 94 ! 95 IF ( MINVAL(champ(:,:)).EQ.MAXVAL(champ(:,:)) .AND. 96 .MINVAL(champ(:,:)).EQ.val_exp ) THEN 97 ! 98 SELECTCASE(varname) 99 ! 100 CASE ('relief') 101 ! 102 ! If we do not have the orography we need to get it 103 ! 104 IF ( .NOT.ALLOCATED(relief)) THEN 105 ! 106 if (present(masque_lu)) then 107 CALL start_init_orog( iml, jml, lon_in, lat_in, 108 . jml2,lon_in2,lat_in2, interbar, masque_lu ) 109 else 110 CALL start_init_orog( iml, jml, lon_in, lat_in, 111 . jml2,lon_in2,lat_in2, interbar) 112 endif 113 ! 114 ENDIF 115 ! 116 IF (SIZE(relief) .NE. SIZE(champ)) THEN 117 ! 118 WRITE(*,*) 'STARTVAR module has been', 119 .' initialized to the wrong size' 120 STOP 121 ! 122 ENDIF 123 ! 124 champ(:,:) = relief(:,:) 125 ! 126 CASE ('rugosite') 127 ! 128 ! If we do not have the orography we need to get it 129 ! 130 IF ( .NOT.ALLOCATED(rugo)) THEN 131 ! 132 CALL start_init_orog( iml, jml, lon_in, lat_in, 133 . jml2,lon_in2,lat_in2 , interbar ) 134 ! 135 ENDIF 136 ! 137 IF (SIZE(rugo) .NE. SIZE(champ)) THEN 138 ! 139 WRITE(*,*) 140 . 'STARTVAR module has been initialized to the wrong size' 141 STOP 142 ! 143 ENDIF 144 ! 145 champ(:,:) = rugo(:,:) 146 ! 147 CASE ('masque') 148 ! 149 ! If we do not have the orography we need to get it 150 ! 151 IF ( .NOT.ALLOCATED(masque)) THEN 152 ! 153 CALL start_init_orog( iml, jml, lon_in, lat_in, 154 . jml2,lon_in2,lat_in2 , interbar ) 155 ! 156 ENDIF 157 ! 158 IF (SIZE(masque) .NE. SIZE(champ)) THEN 159 ! 160 WRITE(*,*) 161 . 'STARTVAR module has been initialized to the wrong size' 162 STOP 163 ! 164 ENDIF 165 ! 166 champ(:,:) = masque(:,:) 167 ! 168 CASE ('surfgeo') 169 ! 170 ! If we do not have the orography we need to get it 171 ! 172 IF ( .NOT.ALLOCATED(phis)) THEN 173 ! 174 CALL start_init_orog( iml, jml, lon_in, lat_in, 175 . jml2,lon_in2, lat_in2 , interbar ) 176 ! 177 ENDIF 178 ! 179 IF (SIZE(phis) .NE. SIZE(champ)) THEN 180 ! 181 WRITE(*,*) 182 . 'STARTVAR module has been initialized to the wrong size' 183 STOP 184 ! 185 ENDIF 186 ! 187 champ(:,:) = phis(:,:) 188 ! 189 CASE ('psol') 190 ! 191 ! If we do not have the orography we need to get it 192 ! 193 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 194 ! 195 CALL start_init_dyn( iml, jml, lon_in, lat_in, 196 . jml2,lon_in2, lat_in2 , interbar ) 197 ! 198 ENDIF 199 ! 200 IF (SIZE(psol_dyn) .NE. SIZE(champ)) THEN 201 ! 202 WRITE(*,*) 203 . 'STARTVAR module has been initialized to the wrong size' 204 STOP 205 ! 206 ENDIF 207 ! 208 champ(:,:) = psol_dyn(:,:) 209 ! 210 CASE DEFAULT 211 ! 212 WRITE(*,*) 'startget_phys2d' 213 WRITE(*,*) 'No rule is present to extract variable', 214 . varname(:LEN_TRIM(varname)),' from any data set' 215 STOP 216 ! 217 END SELECT 218 ! 219 ELSE 220 ! 221 ! There are a few fields we might need if we need to interpolate 3D filed. Thus if they come through here we 222 ! will catch them 223 ! 224 SELECTCASE(varname) 225 ! 226 CASE ('surfgeo') 227 ! 228 IF ( .NOT.ALLOCATED(phis)) THEN 229 ALLOCATE(phis(iml,jml)) 230 ENDIF 231 ! 232 IF (SIZE(phis) .NE. SIZE(champ)) THEN 233 ! 234 WRITE(*,*) 235 . 'STARTVAR module has been initialized to the wrong size' 236 STOP 237 ! 238 ENDIF 239 ! 240 phis(:,:) = champ(:,:) 241 ! 242 END SELECT 243 ! 244 ENDIF 245 ! 246 END SUBROUTINE startget_phys2d 247 ! 248 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 249 ! 250 SUBROUTINE start_init_orog ( iml,jml,lon_in, lat_in,jml2,lon_in2 , 251 , lat_in2 , interbar, masque_lu ) 252 ! 253 INTEGER, INTENT(in) :: iml, jml, jml2 254 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 255 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 256 REAL, intent(in), optional :: masque_lu(iml,jml) 257 LOGICAL interbar 258 ! 259 ! LOCAL 260 ! 261 LOGICAL interbar2 262 REAL :: lev(1), date, dt,chmin,chmax 263 INTEGER :: itau(1), fid 264 INTEGER :: llm_tmp, ttm_tmp 265 INTEGER :: i, j 266 INTEGER :: iret 267 CHARACTER*25 title 268 REAL, ALLOCATABLE :: relief_hi(:,:) 269 REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:) 270 REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:) 271 REAL, ALLOCATABLE :: tmp_var(:,:) 272 INTEGER, ALLOCATABLE :: tmp_int(:,:) 273 ! 274 CHARACTER*120 :: orogfname 275 LOGICAL :: check=.TRUE. 276 ! 277 ! 278 orogfname = 'Relief.nc' 279 ! 280 IF ( check ) WRITE(*,*) 'Reading the high resolution orography' 281 ! 282 CALL flininfo(orogfname,iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) 283 ! 284 ALLOCATE (lat_rel(iml_rel,jml_rel), stat=iret) 285 ALLOCATE (lon_rel(iml_rel,jml_rel), stat=iret) 286 ALLOCATE (relief_hi(iml_rel,jml_rel), stat=iret) 287 ! 288 CALL flinopen(orogfname, .FALSE., iml_rel, jml_rel, 289 .llm_tmp, lon_rel, lat_rel, lev, ttm_tmp, 290 . itau, date, dt, fid) 291 ! 292 CALL flinget(fid, 'RELIEF', iml_rel, jml_rel, llm_tmp, 293 . ttm_tmp, 1, 1, relief_hi) 294 ! 295 CALL flinclo(fid) 296 ! 297 ! In case we have a file which is in degrees we do the transformation 298 ! 299 ALLOCATE(lon_rad(iml_rel)) 300 ALLOCATE(lon_ini(iml_rel)) 301 302 IF ( MAXVAL(lon_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 303 lon_ini(:) = lon_rel(:,1) * 2.0 * ASIN(1.0) / 180.0 304 ELSE 305 lon_ini(:) = lon_rel(:,1) 306 ENDIF 307 308 ALLOCATE(lat_rad(jml_rel)) 309 ALLOCATE(lat_ini(jml_rel)) 310 311 IF ( MAXVAL(lat_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 312 lat_ini(:) = lat_rel(1,:) * 2.0 * ASIN(1.0) / 180.0 313 ELSE 314 lat_ini(:) = lat_rel(1,:) 315 ENDIF 316 ! 317 ! 318 319 title='RELIEF' 320 321 interbar2 = .FALSE. 322 CALL conf_dat2d(title,iml_rel, jml_rel, lon_ini, lat_ini, 323 . lon_rad, lat_rad, relief_hi , interbar2 ) 324 325 IF ( check ) WRITE(*,*) 'Computes all the parameters needed', 326 .' for the gravity wave drag code' 327 ! 328 ! Allocate the data we need to put in the interpolated fields 329 ! 330 ! RELIEF: orographie moyenne 331 ALLOCATE(relief(iml,jml)) 332 ! zphi : orographie moyenne 333 ALLOCATE(phis(iml,jml)) 334 ! zstd: deviation standard de l'orographie sous-maille 335 ALLOCATE(zstd(iml,jml)) 336 ! zsig: pente de l'orographie sous-maille 337 ALLOCATE(zsig(iml,jml)) 338 ! zgam: anisotropy de l'orographie sous maille 339 ALLOCATE(zgam(iml,jml)) 340 ! zthe: orientation de l'axe oriente dans la direction 341 ! de plus grande pente de l'orographie sous maille 342 ALLOCATE(zthe(iml,jml)) 343 ! zpic: hauteur pics de la SSO 344 ALLOCATE(zpic(iml,jml)) 345 ! zval: hauteur vallees de la SSO 346 ALLOCATE(zval(iml,jml)) 347 ! masque : Masque terre ocean 348 ALLOCATE(tmp_int(iml,jml)) 349 ALLOCATE(masque(iml,jml)) 350 351 masque = -99999. 352 if (present(masque_lu)) then 353 masque = masque_lu 354 endif 355 ! 356 CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, 357 . iml-1, jml, lon_in, lat_in, 358 . phis, relief, zstd, zsig, zgam, zthe, zpic, zval, masque) 359 phis = phis * 9.81 360 ! 361 ! masque(:,:) = FLOAT(tmp_int(:,:)) 362 ! 363 ! Compute surface roughness 364 ! 365 IF ( check ) WRITE(*,*) 366 .'Compute surface roughness induced by the orography' 367 ! 368 ALLOCATE(rugo(iml,jml)) 369 ALLOCATE(tmp_var(iml-1,jml)) 370 ! 371 CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, 372 . iml-1, jml, lon_in, lat_in, tmp_var) 373 ! 374 DO j = 1, jml 375 DO i = 1, iml-1 376 rugo(i,j) = tmp_var(i,j) 377 ENDDO 378 rugo(iml,j) = tmp_var(1,j) 379 ENDDO 380 c 381 cc *** rugo n'est pas utilise pour l'instant ****** 382 ! 383 ! Build land-sea mask 384 ! 385 ! 386 RETURN 387 ! 388 END SUBROUTINE start_init_orog 389 ! 390 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 391 ! 392 SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, 393 .lat_in, nbindex, champ, val_exp ,jml2, lon_in2, lat_in2,interbar) 394 ! 395 CHARACTER*(*), INTENT(in) :: varname 396 INTEGER, INTENT(in) :: iml, jml, nbindex, jml2 397 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 398 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 399 REAL, INTENT(inout) :: champ(nbindex) 400 REAL, INTENT(in) :: val_exp 401 LOGICAL interbar 402 ! 403 ! 404 ! This routine only works if the variable does not exist or is constant 405 ! 406 IF ( MINVAL(champ(:)).EQ.MAXVAL(champ(:)) .AND. 407 .MINVAL(champ(:)).EQ.val_exp ) THEN 408 SELECTCASE(varname) 409 CASE ('tsol') 410 IF ( .NOT.ALLOCATED(tsol)) THEN 411 CALL start_init_phys( iml, jml, lon_in, lat_in, 412 . jml2, lon_in2, lat_in2, interbar ) 413 ENDIF 414 IF ( SIZE(tsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 415 WRITE(*,*) 416 . 'STARTVAR module has been initialized to the wrong size' 417 STOP 418 ENDIF 419 CALL gr_dyn_fi(1, iml, jml, nbindex, tsol, champ) 420 CASE ('qsol') 421 IF ( .NOT.ALLOCATED(qsol)) THEN 422 CALL start_init_phys( iml, jml, lon_in, lat_in, 423 . jml2, lon_in2,lat_in2 , interbar ) 424 ENDIF 425 IF ( SIZE(qsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 426 WRITE(*,*) 427 . 'STARTVAR module has been initialized to the wrong size' 428 STOP 429 ENDIF 430 CALL gr_dyn_fi(1, iml, jml, nbindex, qsol, champ) 431 CASE ('psol') 432 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 433 CALL start_init_dyn( iml, jml, lon_in, lat_in, 434 . jml2, lon_in2,lat_in2 , interbar ) 435 ENDIF 436 IF (SIZE(psol_dyn) .NE. SIZE(lon_in)*SIZE(lat_in)) THEN 437 WRITE(*,*) 438 . 'STARTVAR module has been initialized to the wrong size' 439 STOP 440 ENDIF 441 CALL gr_dyn_fi(1, iml, jml, nbindex, psol_dyn, champ) 442 CASE ('zmea') 443 IF ( .NOT.ALLOCATED(relief)) THEN 444 CALL start_init_orog( iml, jml, lon_in, lat_in, 445 . jml2, lon_in2,lat_in2 , interbar ) 446 ENDIF 447 IF ( SIZE(relief) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 448 WRITE(*,*) 449 . 'STARTVAR module has been initialized to the wrong size' 450 STOP 451 ENDIF 452 CALL gr_dyn_fi(1, iml, jml, nbindex, relief, champ) 453 CASE ('zstd') 454 IF ( .NOT.ALLOCATED(zstd)) THEN 455 CALL start_init_orog( iml, jml, lon_in, lat_in, 456 . jml2, lon_in2,lat_in2 , interbar ) 457 ENDIF 458 IF ( SIZE(zstd) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 459 WRITE(*,*) 460 . 'STARTVAR module has been initialized to the wrong size' 461 STOP 462 ENDIF 463 CALL gr_dyn_fi(1, iml, jml, nbindex,zstd, champ) 464 CASE ('zsig') 465 IF ( .NOT.ALLOCATED(zsig)) THEN 466 CALL start_init_orog( iml, jml, lon_in, lat_in, 467 . jml2, lon_in2,lat_in2 , interbar ) 468 ENDIF 469 IF ( SIZE(zsig) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 470 WRITE(*,*) 471 . 'STARTVAR module has been initialized to the wrong size' 472 STOP 473 ENDIF 474 CALL gr_dyn_fi(1, iml, jml, nbindex,zsig, champ) 475 CASE ('zgam') 476 IF ( .NOT.ALLOCATED(zgam)) THEN 477 CALL start_init_orog( iml, jml, lon_in, lat_in, 478 . jml2, lon_in2,lat_in2 , interbar ) 479 ENDIF 480 IF ( SIZE(zgam) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 481 WRITE(*,*) 482 . 'STARTVAR module has been initialized to the wrong size' 483 STOP 484 ENDIF 485 CALL gr_dyn_fi(1, iml, jml, nbindex,zgam, champ) 486 CASE ('zthe') 487 IF ( .NOT.ALLOCATED(zthe)) THEN 488 CALL start_init_orog( iml, jml, lon_in, lat_in, 489 . jml2, lon_in2,lat_in2 , interbar ) 490 ENDIF 491 IF ( SIZE(zthe) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 492 WRITE(*,*) 493 . 'STARTVAR module has been initialized to the wrong size' 494 STOP 495 ENDIF 496 CALL gr_dyn_fi(1, iml, jml, nbindex,zthe, champ) 497 CASE ('zpic') 498 IF ( .NOT.ALLOCATED(zpic)) THEN 499 CALL start_init_orog( iml, jml, lon_in, lat_in, 500 . jml2, lon_in2,lat_in2 , interbar ) 501 ENDIF 502 IF ( SIZE(zpic) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 503 WRITE(*,*) 504 . 'STARTVAR module has been initialized to the wrong size' 505 STOP 506 ENDIF 507 CALL gr_dyn_fi(1, iml, jml, nbindex,zpic, champ) 508 CASE ('zval') 509 IF ( .NOT.ALLOCATED(zval)) THEN 510 CALL start_init_orog( iml, jml, lon_in, lat_in, 511 . jml2, lon_in2,lat_in2 , interbar ) 512 ENDIF 513 IF ( SIZE(zval) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 514 WRITE(*,*) 515 . 'STARTVAR module has been initialized to the wrong size' 516 STOP 517 ENDIF 518 CALL gr_dyn_fi(1, iml, jml, nbindex,zval, champ) 519 CASE ('rads') 520 champ(:) = 0.0 521 CASE ('snow') 522 champ(:) = 0.0 523 cIM "slab" ocean 524 CASE ('tslab') 525 champ(:) = 0.0 526 CASE ('seaice') 527 champ(:) = 0.0 528 CASE ('rugmer') 529 champ(:) = 0.001 530 CASE ('agsno') 531 champ(:) = 50.0 532 CASE DEFAULT 533 WRITE(*,*) 'startget_phys1d' 534 WRITE(*,*) 'No rule is present to extract variable ', 535 . varname(:LEN_TRIM(varname)),' from any data set' 536 STOP 537 END SELECT 538 ELSE 539 ! 540 ! If we see tsol we catch it as we may need it for a 3D interpolation 541 ! 542 SELECTCASE(varname) 543 CASE ('tsol') 544 IF ( .NOT.ALLOCATED(tsol)) THEN 545 ALLOCATE(tsol(SIZE(lon_in),SIZE(lat_in) )) 546 ENDIF 547 CALL gr_fi_dyn(1, iml, jml, nbindex, champ, tsol) 52 USE ioipsl 53 IMPLICIT NONE 54 55 PRIVATE 56 PUBLIC startget 57 INTERFACE startget 58 MODULE PROCEDURE startget_phys1d, startget_phys2d, startget_dyn 59 END INTERFACE 60 61 REAL, SAVE :: deg2rad, pi 62 INTEGER, SAVE :: iml_rel, jml_rel 63 INTEGER, SAVE :: fid_phys, iml_phys, jml_phys 64 INTEGER, SAVE :: fid_dyn, iml_dyn, jml_dyn, llm_dyn, ttm_dyn 65 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lon_phys, lon_dyn 66 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lat_phys, lat_dyn 67 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lon_rug, lon_alb, lon_rel 68 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lat_rug, lat_alb, lat_rel 69 REAL, DIMENSION(:), ALLOCATABLE, TARGET, SAVE :: levdyn_ini 70 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: relief, zstd, zsig, zgam 71 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: masque, zthe, zpic, zval 72 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: rugo, phis, tsol, qsol 73 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: psol_dyn 74 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET, SAVE :: var_ana3d 75 76 CONTAINS 77 78 !------------------------------------------------------------------------------- 79 ! 80 SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, lat_in, nbindex, champ, & 81 val_exp ,jml2, lon_in2, lat_in2, ibar) 82 ! 83 !------------------------------------------------------------------------------- 84 ! Comment: 85 ! This routine only works if the variable does not exist or is constant. 86 !------------------------------------------------------------------------------- 87 ! Arguments: 88 CHARACTER(LEN=*), INTENT(IN) :: varname 89 INTEGER, INTENT(IN) :: iml, jml 90 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 91 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 92 INTEGER, INTENT(IN) :: nbindex 93 REAL, DIMENSION(nbindex), INTENT(INOUT) :: champ 94 REAL, INTENT(IN) :: val_exp 95 INTEGER, INTENT(IN) :: jml2 96 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 97 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 98 LOGICAL, INTENT(IN) :: ibar 99 !------------------------------------------------------------------------------- 100 ! Local variables: 101 #include "iniprint.h" 102 REAL, DIMENSION(:,:), POINTER :: v2d 103 !------------------------------------------------------------------------------- 104 v2d=>NULL() 105 IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN 106 107 !--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES 108 SELECT CASE(varname) 109 CASE('tsol') 110 IF(.NOT.ALLOCATED(tsol)) & 111 CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 112 CASE('qsol') 113 IF(.NOT.ALLOCATED(qsol)) & 114 CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 115 CASE('psol') 116 IF(.NOT.ALLOCATED(psol_dyn)) & 117 CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 118 CASE('zmea','zstd','zsig','zgam','zthe','zpic','zval') 119 IF(.NOT.ALLOCATED(relief)) & 120 CALL start_init_orog(iml,jml,lon_in,lat_in) 121 CASE('rads','snow','tslab','seaice','rugmer','agsno') 122 CASE DEFAULT 123 WRITE(lunout,*)'startget_phys1d' 124 WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)& 125 //' from any data set'; STOP 126 END SELECT 127 128 !--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE 129 SELECT CASE(varname) 130 CASE('rads','snow','tslab','seaice'); champ=0.0 131 CASE('rugmer'); champ(:)=0.001 132 CASE('agsno'); champ(:)=50.0 133 CASE DEFAULT 134 SELECT CASE(varname) 135 CASE('tsol'); v2d=>tsol 136 CASE('qsol'); v2d=>qsol 137 CASE('psol'); v2d=>psol_dyn 138 CASE('zmea'); v2d=>relief 139 CASE('zstd'); v2d=>zstd 140 CASE('zsig'); v2d=>zsig 141 CASE('zgam'); v2d=>zgam 142 CASE('zthe'); v2d=>zthe 143 CASE('zpic'); v2d=>zpic 144 CASE('zval'); v2d=>zval 548 145 END SELECT 549 ENDIF 550 END SUBROUTINE startget_phys1d 551 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 552 ! 553 SUBROUTINE start_init_phys( iml, jml, lon_in, lat_in, jml2, 554 . lon_in2, lat_in2 , interbar ) 555 ! 556 INTEGER, INTENT(in) :: iml, jml ,jml2 557 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 558 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 559 LOGICAL interbar 560 ! 561 ! LOCAL 562 ! 563 !ac REAL :: lev(1), date, dt 564 REAL :: date, dt 565 REAL, DIMENSION(:), ALLOCATABLE :: levphys_ini 566 !ac 567 INTEGER :: itau(1) 568 INTEGER :: llm_tmp, ttm_tmp 569 INTEGER :: i, j 570 ! 571 CHARACTER*25 title 572 CHARACTER*120 :: physfname 573 LOGICAL :: check=.TRUE. 574 ! 575 REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:) 576 REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:) 577 REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:) 578 ! 579 physfname = 'ECPHY.nc' 580 ! 581 IF ( check ) WRITE(*,*) 'Opening the surface analysis' 582 ! 583 CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp, 584 . ttm_tmp, fid_phys) 585 ! 586 ALLOCATE (lat_phys(iml_phys,jml_phys)) 587 ALLOCATE (lon_phys(iml_phys,jml_phys)) 588 !ac 589 ALLOCATE (levphys_ini(llm_tmp)) 590 ! 591 ! CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, 592 ! . llm_tmp, lon_phys, lat_phys, lev, ttm_tmp, 593 ! . itau, date, dt, fid_phys) 594 ! 595 CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, 596 . llm_tmp, lon_phys, lat_phys, levphys_ini, ttm_tmp, 597 . itau, date, dt, fid_phys) 598 ! 599 DEALLOCATE (levphys_ini) 600 !ac 601 ! 602 ! Allocate the space we will need to get the data out of this file 603 ! 604 ALLOCATE(var_ana(iml_phys, jml_phys)) 605 ! 606 ! In case we have a file which is in degrees we do the transformation 607 ! 608 ALLOCATE(lon_rad(iml_phys)) 609 ALLOCATE(lon_ini(iml_phys)) 610 611 IF ( MAXVAL(lon_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 612 lon_ini(:) = lon_phys(:,1) * 2.0 * ASIN(1.0) / 180.0 613 ELSE 614 lon_ini(:) = lon_phys(:,1) 615 ENDIF 616 617 ALLOCATE(lat_rad(jml_phys)) 618 ALLOCATE(lat_ini(jml_phys)) 619 620 IF ( MAXVAL(lat_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 621 lat_ini(:) = lat_phys(1,:) * 2.0 * ASIN(1.0) / 180.0 622 ELSE 623 lat_ini(:) = lat_phys(1,:) 624 ENDIF 625 626 627 ! 628 ! We get the two standard varibales 629 ! Surface temperature 630 ! 631 ALLOCATE(tsol(iml,jml)) 632 ALLOCATE(tmp_var(iml-1,jml)) 633 ! 634 ! 635 636 CALL flinget(fid_phys, 'ST', iml_phys, jml_phys, 637 .llm_tmp, ttm_tmp, 1, 1, var_ana) 638 639 title='ST' 640 CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, 641 . lon_rad, lat_rad, var_ana , interbar ) 642 643 IF ( interbar ) THEN 644 WRITE(6,*) '-------------------------------------------------', 645 ,'--------------' 646 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 647 , ' pour ST $$$ ' 648 WRITE(6,*) '-------------------------------------------------', 649 ,'--------------' 650 CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad , 651 , var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var ) 652 ELSE 653 CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad, 654 . var_ana, iml-1, jml, lon_in, lat_in, tmp_var ) 655 ENDIF 656 657 CALL gr_int_dyn(tmp_var, tsol, iml-1, jml) 658 ! 659 ! Soil moisture 660 ! 661 ALLOCATE(qsol(iml,jml)) 662 CALL flinget(fid_phys, 'CDSW', iml_phys, jml_phys, 663 . llm_tmp, ttm_tmp, 1, 1, var_ana) 664 665 title='CDSW' 666 CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, 667 . lon_rad, lat_rad, var_ana, interbar ) 668 669 IF ( interbar ) THEN 670 WRITE(6,*) '-------------------------------------------------', 671 ,'--------------' 672 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 673 , ' pour CDSW $$$ ' 674 WRITE(6,*) '-------------------------------------------------', 675 ,'--------------' 676 CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad , 677 , var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var ) 678 ELSE 679 CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad, 680 . var_ana, iml-1, jml, lon_in, lat_in, tmp_var ) 681 ENDIF 682 c 683 CALL gr_int_dyn(tmp_var, qsol, iml-1, jml) 684 ! 685 CALL flinclo(fid_phys) 686 ! 687 END SUBROUTINE start_init_phys 688 ! 689 ! 690 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 691 ! 692 ! 693 SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in, 694 . lml, pls, workvar, champ, val_exp,jml2, lon_in2, lat_in2 , 695 , interbar ) 696 ! 697 ! ARGUMENTS 698 ! 699 CHARACTER*(*), INTENT(in) :: varname 700 INTEGER, INTENT(in) :: iml, jml, lml, jml2 701 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 702 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 703 REAL, INTENT(in) :: pls(iml, jml, lml) 704 REAL, INTENT(in) :: workvar(iml, jml, lml) 705 REAL, INTENT(inout) :: champ(iml, jml, lml) 706 REAL, INTENT(in) :: val_exp 707 LOGICAL interbar 708 ! 709 ! LOCAL 710 ! 711 INTEGER :: il, ij, ii 712 REAL :: xppn, xpps 713 ! 146 IF(SIZE(v2d)/=SIZE(lon_in)*SIZE(lat_in)) THEN 147 WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size' 148 STOP 149 END IF 150 CALL gr_dyn_fi(1,iml,jml,nbindex,v2d,champ) 151 END SELECT 152 153 ELSE 154 155 !--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION 156 SELECT CASE(varname) 157 CASE('tsol') 158 IF(.NOT.ALLOCATED(tsol)) ALLOCATE(tsol(iml,jml)) 159 CALL gr_fi_dyn(1,iml,jml,nbindex,champ,tsol) 160 END SELECT 161 162 END IF 163 164 END SUBROUTINE startget_phys1d 165 ! 166 !------------------------------------------------------------------------------- 167 168 169 !------------------------------------------------------------------------------- 170 ! 171 SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in, champ, val_exp, & 172 jml2, lon_in2, lat_in2 , ibar, msk) 173 ! 174 !------------------------------------------------------------------------------- 175 ! Comment: 176 ! This routine only works if the variable does not exist or is constant. 177 !------------------------------------------------------------------------------- 178 ! Arguments: 179 CHARACTER(LEN=*), INTENT(IN) :: varname 180 INTEGER, INTENT(IN) :: iml, jml 181 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 182 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 183 REAL, DIMENSION(iml,jml), INTENT(INOUT) :: champ 184 REAL, INTENT(IN) :: val_exp 185 INTEGER, INTENT(IN) :: jml2 186 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 187 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 188 LOGICAL, INTENT(IN) :: ibar 189 REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: msk 190 !------------------------------------------------------------------------------- 191 ! Local variables: 192 #include "iniprint.h" 193 REAL, DIMENSION(:,:), POINTER :: v2d=>NULL() 194 LOGICAL :: lrelief1, lrelief2 195 !------------------------------------------------------------------------------- 196 v2d=>NULL() 197 lrelief1=(.NOT.ALLOCATED(relief).AND. PRESENT(msk)) 198 lrelief2=(.NOT.ALLOCATED(relief).AND..NOT.PRESENT(msk)) 199 IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN 200 201 !--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES 202 SELECT CASE(varname) 203 CASE('psol') 204 IF(.NOT.ALLOCATED(psol_dyn)) & 205 CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 206 CASE('relief') 207 IF(lrelief1) CALL start_init_orog(iml,jml,lon_in,lat_in,msk) 208 IF(lrelief2) CALL start_init_orog(iml,jml,lon_in,lat_in) 209 CASE('surfgeo') 210 IF(.NOT.ALLOCATED(phis)) CALL start_init_orog(iml,jml,lon_in,lat_in) 211 CASE('rugosite','masque') 212 IF(.NOT.ALLOCATED(rugo)) CALL start_init_orog(iml,jml,lon_in,lat_in) 213 CASE DEFAULT 214 WRITE(lunout,*)'startget_phys2d' 215 WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)& 216 //' from any data set'; STOP 217 END SELECT 218 219 !--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE 220 SELECT CASE(varname) 221 CASE('psol'); v2d=>psol_dyn 222 CASE('relief'); v2d=>relief 223 CASE('rugosite'); v2d=>rugo 224 CASE('masque'); v2d=>masque 225 CASE('surfgeo'); v2d=>phis 226 END SELECT 227 IF(SIZE(champ)/=SIZE(v2d)) THEN 228 WRITE(lunout,*) 'STARTVAR module has been initialized to the wrong size' 229 STOP 230 END IF 231 232 champ(:,:)=v2d(:,:) 233 234 ELSE 235 236 !--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION 237 SELECT CASE(varname) 238 CASE ('surfgeo') 239 IF(.NOT.ALLOCATED(phis)) ALLOCATE(phis(iml,jml)) 240 IF(SIZE(phis)/=SIZE(champ)) THEN 241 WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size' 242 STOP 243 END IF 244 phis(:,:)=champ(:,:) 245 END SELECT 246 247 END IF 248 249 END SUBROUTINE startget_phys2d 250 ! 251 !------------------------------------------------------------------------------- 252 253 254 !------------------------------------------------------------------------------- 255 ! 256 SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in, lml, pls,workvar,& 257 champ, val_exp, jml2, lon_in2, lat_in2, ibar) 258 !------------------------------------------------------------------------------- 259 ! Comment: 260 ! This routine only works if the variable does not exist or is constant. 261 !------------------------------------------------------------------------------- 262 ! Arguments: 263 CHARACTER(LEN=*), INTENT(IN) :: varname 264 INTEGER, INTENT(IN) :: iml, jml 265 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 266 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 267 INTEGER, INTENT(IN) :: lml 268 REAL, DIMENSION(iml,jml,lml), INTENT(IN) :: pls, workvar 269 REAL, DIMENSION(iml,jml,lml), INTENT(INOUT) :: champ 270 REAL, INTENT(IN) :: val_exp 271 INTEGER, INTENT(IN) :: jml2 272 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 273 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 274 LOGICAL, INTENT(IN) :: ibar 275 !------------------------------------------------------------------------------- 276 ! Local variables: 277 #include "iniprint.h" 278 #include "dimensions.h" 279 #include "comconst.h" 280 #include "paramet.h" 281 #include "comgeom2.h" 282 REAL, DIMENSION(:,:,:), POINTER :: v3d=>NULL() 283 CHARACTER(LEN=10) :: vname 284 INTEGER :: il 285 REAL :: xppn, xpps 286 !------------------------------------------------------------------------------- 287 NULLIFY(v3d) 288 IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN 289 290 !--- READING UNALLOCATED FILES 291 IF(.NOT.ALLOCATED(psol_dyn)) & 292 CALL start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 293 294 !--- CHECKING IF THE FIELD IS KNOWN AND INTERPOLATING 3D FIELDS 295 SELECT CASE(varname) 296 CASE('u'); vname='U' 297 CASE('v'); vname='V' 298 CASE('t','tpot'); vname='TEMP' 299 CASE('q'); vname='R' 300 CASE DEFAULT 301 WRITE(lunout,*)'startget_dyn' 302 WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)& 303 //' from any data set'; STOP 304 END SELECT 305 CALL start_inter_3d(TRIM(vname), iml, jml, lml, lon_in, lat_in, jml2, & 306 lon_in2, lat_in2, pls, champ,ibar ) 307 308 !--- COMPUTING THE REQUIRED FILED 309 SELECT CASE(varname) 310 CASE('u') !--- Eastward wind 311 DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO 312 champ(iml,:,:)=champ(1,:,:) 313 314 CASE('v') !--- Northward wind 315 DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO 316 champ(iml,:,:)=champ(1,:,:) 317 318 CASE('tpot') !--- Temperature 319 IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN 320 champ=champ*cpp/workvar 321 DO il=1,lml 322 xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln 323 xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols 324 champ(:,1 ,il) = xppn 325 champ(:,jml,il) = xpps 326 END DO 327 ELSE 328 WRITE(lunout,*)'Could not compute potential temperature as the' 329 WRITE(lunout,*)'Exner function is missing or constant.'; STOP 330 END IF 331 332 CASE('q') !--- Relat. humidity 333 IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN 334 champ=0.01*champ*workvar 335 WHERE(champ<0.) champ=1.0E-10 336 DO il=1,lml 337 xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln 338 xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols 339 champ(:,1 ,il) = xppn 340 champ(:,jml,il) = xpps 341 END DO 342 ELSE 343 WRITE(lunout,*)'Could not compute specific humidity as the' 344 WRITE(lunout,*)'saturated humidity is missing or constant.'; STOP 345 END IF 346 347 END SELECT 348 349 END IF 350 351 END SUBROUTINE startget_dyn 352 ! 353 !------------------------------------------------------------------------------- 354 355 356 !------------------------------------------------------------------------------- 357 ! 358 SUBROUTINE start_init_orog(iml,jml,lon_in,lat_in,masque_lu) 359 ! 360 !------------------------------------------------------------------------------- 361 ! Arguments: 362 INTEGER, INTENT(IN) :: iml, jml 363 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 364 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 365 REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: masque_lu 366 !------------------------------------------------------------------------------- 367 ! Local variables: 368 #include "iniprint.h" 369 CHARACTER(LEN=25) :: title 370 CHARACTER(LEN=120) :: orofname 371 LOGICAL :: check=.TRUE. 372 REAL, DIMENSION(1) :: lev 373 REAL :: date, dt 374 INTEGER, DIMENSION(1) :: itau 375 INTEGER :: fid, llm_tmp, ttm_tmp 376 REAL, DIMENSION(:,:), ALLOCATABLE :: relief_hi, tmp_var 377 REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini 378 !------------------------------------------------------------------------------- 379 pi=2.0*ASIN(1.0); deg2rad=pi/180.0 380 381 orofname = 'Relief.nc'; title='RELIEF' 382 IF(check) WRITE(lunout,*)'Reading the high resolution orography' 383 CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) 384 385 ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel)) 386 CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,& 387 lev, ttm_tmp, itau, date, dt, fid) 388 ALLOCATE(relief_hi(iml_rel,jml_rel)) 389 CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi) 390 CALL flinclo(fid) 391 392 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 393 ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel)) 394 lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad 395 lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad 396 397 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS 398 ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel)) 399 CALL conf_dat2d(title, iml_rel, jml_rel, lon_ini, lat_ini, lon_rad, lat_rad, & 400 relief_hi, .FALSE.) 401 DEALLOCATE(lon_ini,lat_ini) 402 403 !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro 404 IF(check) WRITE(lunout,*)'Computes all parameters needed for gravity wave dra& 405 &g code' 406 407 ALLOCATE(phis(iml,jml)) ! Geopotentiel au sol 408 ALLOCATE(zstd(iml,jml)) ! Deviation standard de l'orographie sous-maille 409 ALLOCATE(zsig(iml,jml)) ! Pente de l'orographie sous-maille 410 ALLOCATE(zgam(iml,jml)) ! Anisotropie de l'orographie sous maille 411 ALLOCATE(zthe(iml,jml)) ! Orientation axe +grande pente d'oro sous maille 412 ALLOCATE(zpic(iml,jml)) ! Hauteur pics de la SSO 413 ALLOCATE(zval(iml,jml)) ! Hauteur vallees de la SSO 414 ALLOCATE(relief(iml,jml)) ! Orographie moyenne 415 ALLOCATE(masque(iml,jml)) ! Masque terre ocean 416 masque = -99999. 417 IF(PRESENT(masque_lu)) masque=masque_lu 418 419 CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml, & 420 lon_in, lat_in, phis, relief, zstd, zsig, zgam, zthe, zpic, zval, masque) 421 phis = phis * 9.81 422 423 !--- SURFACE ROUGHNESS COMPUTATION (UNUSED FOR THE MOMENT !!! ) 424 IF(check) WRITE(lunout,*)'Compute surface roughness induced by the orography' 425 ALLOCATE(rugo (iml ,jml)) 426 ALLOCATE(tmp_var(iml-1,jml)) 427 CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml, & 428 lon_in, lat_in, tmp_var) 429 rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:) 430 DEALLOCATE(relief_hi,tmp_var,lon_rad,lat_rad) 431 RETURN 432 433 END SUBROUTINE start_init_orog 434 ! 435 !------------------------------------------------------------------------------- 436 437 438 !------------------------------------------------------------------------------- 439 ! 440 SUBROUTINE start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 441 ! 442 !------------------------------------------------------------------------------- 443 ! Arguments: 444 INTEGER, INTENT(IN) :: iml, jml 445 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 446 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 447 INTEGER, INTENT(IN) :: jml2 448 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 449 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 450 LOGICAL, INTENT(IN) :: ibar 451 !------------------------------------------------------------------------------- 452 ! Local variables: 453 #include "iniprint.h" 454 CHARACTER(LEN=25) :: title 455 CHARACTER(LEN=120) :: physfname 456 LOGICAL :: check=.TRUE. 457 REAL :: date, dt 458 INTEGER, DIMENSION(1) :: itau 459 INTEGER :: llm_tmp, ttm_tmp 460 REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana 461 REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini 462 REAL, DIMENSION(:), ALLOCATABLE :: levphys_ini 463 !------------------------------------------------------------------------------- 464 physfname = 'ECPHY.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0 465 IF(check) WRITE(lunout,*)'Opening the surface analysis' 466 CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp, ttm_tmp, fid_phys) 467 468 ALLOCATE(lat_phys(iml_phys,jml_phys)) 469 ALLOCATE(lon_phys(iml_phys,jml_phys)) 470 ALLOCATE(levphys_ini(llm_tmp)) 471 CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_tmp, lon_phys, & 472 lat_phys, levphys_ini, ttm_tmp, itau, date, dt, fid_phys) 473 DEALLOCATE(levphys_ini) 474 475 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 476 ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys)) 477 lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad 478 lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad 479 480 ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys)) 481 482 !--- SURFACE TEMPERATURE 483 title='ST' 484 ALLOCATE(tsol(iml,jml)) 485 CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana) 486 CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,& 487 var_ana , ibar ) 488 CALL interp_startvar(title, ibar, .TRUE., & 489 iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & 490 lon_in, lat_in, lon_in2, lat_in2, tsol) 491 492 !--- SOIL MOISTURE 493 title='CDSW' 494 ALLOCATE(qsol(iml,jml)) 495 CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana) 496 CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,& 497 var_ana, ibar ) 498 CALL interp_startvar(title, ibar, .TRUE., & 499 iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & 500 lon_in, lat_in, lon_in2, lat_in2, qsol) 501 502 CALL flinclo(fid_phys) 503 504 DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) 505 506 END SUBROUTINE start_init_phys 507 ! 508 !------------------------------------------------------------------------------- 509 510 511 !------------------------------------------------------------------------------- 512 ! 513 SUBROUTINE start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 514 ! 515 !------------------------------------------------------------------------------- 516 ! Arguments: 517 INTEGER, INTENT(IN) :: iml, jml 518 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 519 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 520 INTEGER, INTENT(IN) :: jml2 521 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 522 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 523 LOGICAL, INTENT(IN) :: ibar 524 !------------------------------------------------------------------------------- 525 ! Local variables: 526 #include "iniprint.h" 714 527 #include "dimensions.h" 715 528 #include "paramet.h" 716 529 #include "comgeom2.h" 717 #include "comconst.h" 718 ! 719 ! This routine only works if the variable does not exist or is constant 720 ! 721 IF ( MINVAL(champ(:,:,:)).EQ.MAXVAL(champ(:,:,:)) .AND. 722 . MINVAL(champ(:,:,:)).EQ.val_exp ) THEN 723 ! 724 SELECTCASE(varname) 725 CASE ('u') 726 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 727 CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 , 728 . lon_in2,lat_in2 , interbar ) 729 ENDIF 730 CALL start_inter_3d('U', iml, jml, lml, lon_in, 731 . lat_in, jml2, lon_in2, lat_in2, pls, champ,interbar ) 732 DO il=1,lml 733 DO ij=1,jml 734 DO ii=1,iml-1 735 champ(ii,ij,il) = champ(ii,ij,il) * cu(ii,ij) 736 ENDDO 737 champ(iml,ij, il) = champ(1,ij, il) 738 ENDDO 739 ENDDO 740 CASE ('v') 741 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 742 CALL start_init_dyn( iml, jml, lon_in, lat_in , jml2, 743 . lon_in2, lat_in2 , interbar ) 744 ENDIF 745 CALL start_inter_3d('V', iml, jml, lml, lon_in, 746 . lat_in, jml2, lon_in2, lat_in2, pls, champ, interbar ) 747 DO il=1,lml 748 DO ij=1,jml 749 DO ii=1,iml-1 750 champ(ii,ij,il) = champ(ii,ij,il) * cv(ii,ij) 751 ENDDO 752 champ(iml,ij, il) = champ(1,ij, il) 753 ENDDO 754 ENDDO 755 CASE ('t') 756 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 757 CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 , 758 . lon_in2, lat_in2 ,interbar ) 759 ENDIF 760 CALL start_inter_3d('TEMP', iml, jml, lml, lon_in, 761 . lat_in, jml2, lon_in2, lat_in2, pls, champ, interbar ) 762 763 CASE ('tpot') 764 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 765 CALL start_init_dyn( iml, jml, lon_in, lat_in , jml2 , 766 . lon_in2, lat_in2 , interbar ) 767 ENDIF 768 CALL start_inter_3d('TEMP', iml, jml, lml, lon_in, 769 . lat_in, jml2, lon_in2, lat_in2, pls, champ, interbar ) 770 IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) ) 771 . THEN 772 DO il=1,lml 773 DO ij=1,jml 774 DO ii=1,iml-1 775 champ(ii,ij,il) = champ(ii,ij,il) * cpp 776 . / workvar(ii,ij,il) 777 ENDDO 778 champ(iml,ij,il) = champ(1,ij,il) 779 ENDDO 780 ENDDO 781 DO il=1,lml 782 xppn = SUM(aire(:,1)*champ(:,1,il))/apoln 783 xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols 784 champ(:,1,il) = xppn 785 champ(:,jml,il) = xpps 786 ENDDO 787 ELSE 788 WRITE(*,*)'Could not compute potential temperature as the' 789 WRITE(*,*)'Exner function is missing or constant.' 790 STOP 791 ENDIF 792 CASE ('q') 793 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 794 CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 , 795 . lon_in2, lat_in2 , interbar ) 796 ENDIF 797 CALL start_inter_3d('R', iml, jml, lml, lon_in, lat_in, 798 . jml2, lon_in2, lat_in2, pls, champ, interbar ) 799 IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) ) 800 . THEN 801 DO il=1,lml 802 DO ij=1,jml 803 DO ii=1,iml-1 804 champ(ii,ij,il) = 0.01 * champ(ii,ij,il) * 805 . workvar(ii,ij,il) 806 ENDDO 807 champ(iml,ij,il) = champ(1,ij,il) 808 ENDDO 809 ENDDO 810 WHERE ( champ .LT. 0.) champ = 1.0E-10 811 DO il=1,lml 812 xppn = SUM(aire(:,1)*champ(:,1,il))/apoln 813 xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols 814 champ(:,1,il) = xppn 815 champ(:,jml,il) = xpps 816 ENDDO 817 ELSE 818 WRITE(*,*)'Could not compute specific humidity as the' 819 WRITE(*,*)'saturated humidity is missing or constant.' 820 STOP 821 ENDIF 822 CASE DEFAULT 823 WRITE(*,*) 'startget_dyn' 824 WRITE(*,*) 'No rule is present to extract variable ', 825 . varname(:LEN_TRIM(varname)),' from any data set' 826 STOP 827 END SELECT 828 ENDIF 829 END SUBROUTINE startget_dyn 830 ! 831 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 832 ! 833 SUBROUTINE start_init_dyn( iml, jml, lon_in, lat_in,jml2,lon_in2 , 834 , lat_in2 , interbar ) 835 ! 836 INTEGER, INTENT(in) :: iml, jml, jml2 837 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 838 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 839 LOGICAL interbar 840 ! 841 ! LOCAL 842 ! 843 REAL :: lev(1), date, dt 844 INTEGER :: itau(1) 845 INTEGER :: i, j 846 integer :: iret 847 ! 848 CHARACTER*120 :: physfname 849 LOGICAL :: check=.TRUE. 850 ! 851 REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:) 852 REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:) 853 REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:), z(:,:) 854 REAL, ALLOCATABLE :: xppn(:), xpps(:) 855 LOGICAL :: allo 856 ! 857 ! 858 #include "dimensions.h" 859 #include "paramet.h" 860 #include "comgeom2.h" 861 862 CHARACTER*25 title 863 864 ! 865 physfname = 'ECDYN.nc' 866 ! 867 IF ( check ) WRITE(*,*) 'Opening the surface analysis' 868 ! 869 CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, 870 . ttm_dyn, fid_dyn) 871 IF ( check ) WRITE(*,*) 'Values read: ', iml_dyn, jml_dyn, 872 . llm_dyn, ttm_dyn 873 ! 874 ALLOCATE (lat_dyn(iml_dyn,jml_dyn), stat=iret) 875 ALLOCATE (lon_dyn(iml_dyn,jml_dyn), stat=iret) 876 ALLOCATE (levdyn_ini(llm_dyn), stat=iret) 877 ! 878 CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, 879 . lon_dyn, lat_dyn, levdyn_ini, ttm_dyn, 880 . itau, date, dt, fid_dyn) 881 ! 882 883 allo = allocated (var_ana) 884 if (allo) then 885 DEALLOCATE(var_ana, stat=iret) 886 endif 887 ALLOCATE(var_ana(iml_dyn, jml_dyn), stat=iret) 888 889 allo = allocated (lon_rad) 890 if (allo) then 891 DEALLOCATE(lon_rad, stat=iret) 892 endif 893 894 ALLOCATE(lon_rad(iml_dyn), stat=iret) 895 ALLOCATE(lon_ini(iml_dyn)) 896 897 IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 898 lon_ini(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0 899 ELSE 900 lon_ini(:) = lon_dyn(:,1) 901 ENDIF 902 903 ALLOCATE(lat_rad(jml_dyn)) 904 ALLOCATE(lat_ini(jml_dyn)) 905 906 IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 907 lat_ini(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0 908 ELSE 909 lat_ini(:) = lat_dyn(1,:) 910 ENDIF 911 ! 912 913 914 ALLOCATE(z(iml, jml)) 915 ALLOCATE(tmp_var(iml-1,jml)) 916 ! 917 CALL flinget(fid_dyn, 'Z', iml_dyn, jml_dyn, 0, ttm_dyn, 918 . 1, 1, var_ana) 919 c 920 title='Z' 921 CALL conf_dat2d( title,iml_dyn, jml_dyn,lon_ini, lat_ini, 922 . lon_rad, lat_rad, var_ana, interbar ) 923 c 924 IF ( interbar ) THEN 925 WRITE(6,*) '-------------------------------------------------', 926 ,'--------------' 927 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 928 , ' pour Z $$$ ' 929 WRITE(6,*) '-------------------------------------------------', 930 ,'--------------' 931 CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad , 932 , var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var) 933 ELSE 934 CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana, 935 . iml-1, jml, lon_in, lat_in, tmp_var) 936 ENDIF 937 938 CALL gr_int_dyn(tmp_var, z, iml-1, jml) 939 ! 940 ALLOCATE(psol_dyn(iml, jml)) 941 ! 942 CALL flinget(fid_dyn, 'SP', iml_dyn, jml_dyn, 0, ttm_dyn, 943 . 1, 1, var_ana) 944 945 title='SP' 946 CALL conf_dat2d( title,iml_dyn, jml_dyn,lon_ini, lat_ini, 947 . lon_rad, lat_rad, var_ana, interbar ) 948 949 IF ( interbar ) THEN 950 WRITE(6,*) '-------------------------------------------------', 951 ,'--------------' 952 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 953 , ' pour SP $$$ ' 954 WRITE(6,*) '-------------------------------------------------', 955 ,'--------------' 956 CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad , 957 , var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var) 958 ELSE 959 CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana, 960 . iml-1, jml, lon_in, lat_in, tmp_var ) 961 ENDIF 962 963 CALL gr_int_dyn(tmp_var, psol_dyn, iml-1, jml) 964 ! 965 IF ( .NOT.ALLOCATED(tsol)) THEN 966 ! These variables may have been allocated by the need to 967 ! create a start field for them or by the varibale 968 ! coming out of the restart file. In case we dor have it we will initialize it. 969 ! 970 CALL start_init_phys( iml, jml, lon_in, lat_in,jml2,lon_in2, 971 . lat_in2 , interbar ) 972 ELSE 973 IF ( SIZE(tsol) .NE. SIZE(psol_dyn) ) THEN 974 WRITE(*,*) 'start_init_dyn :' 975 WRITE(*,*) 'The temperature field we have does not ', 976 . 'have the right size' 977 STOP 978 ENDIF 979 ENDIF 980 IF ( .NOT.ALLOCATED(phis)) THEN 981 ! 982 ! These variables may have been allocated by the need to create a start field for them or by the varibale 983 ! coming out of the restart file. In case we dor have it we will initialize it. 984 ! 985 CALL start_init_orog( iml, jml, lon_in, lat_in, jml2, lon_in2 , 986 . lat_in2 , interbar ) 987 ! 988 ELSE 989 ! 990 IF (SIZE(phis) .NE. SIZE(psol_dyn)) THEN 991 ! 992 WRITE(*,*) 'start_init_dyn :' 993 WRITE(*,*) 'The orography field we have does not ', 994 . ' have the right size' 995 STOP 996 ENDIF 997 ! 998 ENDIF 999 ! 1000 ! PSOL is computed in Pascals 1001 ! 1002 ! 1003 DO j = 1, jml 1004 DO i = 1, iml-1 1005 psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j)) 1006 . /287.0/tsol(i,j)) 1007 ENDDO 1008 psol_dyn(iml,j) = psol_dyn(1,j) 1009 ENDDO 1010 ! 1011 ! 1012 ALLOCATE(xppn(iml-1)) 1013 ALLOCATE(xpps(iml-1)) 1014 ! 1015 DO i = 1, iml-1 1016 xppn(i) = aire( i,1) * psol_dyn( i,1) 1017 xpps(i) = aire( i,jml) * psol_dyn( i,jml) 1018 ENDDO 1019 ! 1020 DO i = 1, iml 1021 psol_dyn(i,1 ) = SUM(xppn)/apoln 1022 psol_dyn(i,jml) = SUM(xpps)/apols 1023 ENDDO 1024 ! 1025 RETURN 1026 ! 1027 END SUBROUTINE start_init_dyn 1028 ! 1029 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1030 ! 1031 SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, 1032 . lat_in, jml2, lon_in2, lat_in2, pls_in, var3d, interbar ) 1033 ! 1034 ! This subroutine gets a variables from a 3D file and does the interpolations needed 1035 ! 1036 ! 1037 ! ARGUMENTS 1038 ! 1039 CHARACTER*(*) :: varname 1040 INTEGER :: iml, jml, lml, jml2 1041 REAL :: lon_in(iml), lat_in(jml), pls_in(iml, jml, lml) 1042 REAL :: lon_in2(iml) , lat_in2(jml2) 1043 REAL :: var3d(iml, jml, lml) 1044 LOGICAL interbar 1045 real chmin,chmax 1046 ! 1047 ! LOCAL 1048 ! 1049 CHARACTER*25 title 1050 INTEGER :: ii, ij, il, jsort,i,j,l 1051 REAL :: bx, by 1052 REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:) 1053 REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:) , lev_dyn(:) 1054 REAL, ALLOCATABLE :: var_tmp2d(:,:), var_tmp3d(:,:,:) 1055 REAL, ALLOCATABLE :: ax(:), ay(:), yder(:) 1056 ! REAL, ALLOCATABLE :: varrr(:,:,:) 1057 INTEGER, ALLOCATABLE :: lind(:) 1058 ! 1059 LOGICAL :: check = .TRUE. 1060 ! 1061 IF ( .NOT. ALLOCATED(var_ana3d)) THEN 1062 ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn)) 1063 ENDIF 1064 ! ALLOCATE(varrr(iml_dyn, jml_dyn, llm_dyn)) 1065 ! 1066 ! 1067 IF ( check) WRITE(*,*) 'Going into flinget to extract the 3D ', 1068 . ' field.', fid_dyn 1069 IF ( check) WRITE(*,*) fid_dyn, varname, iml_dyn, jml_dyn, 1070 . llm_dyn,ttm_dyn 1071 ! 1072 CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, 1073 . ttm_dyn, 1, 1, var_ana3d) 1074 ! 1075 IF ( check) WRITE(*,*) 'Allocating space for the interpolation', 1076 . iml, jml, llm_dyn 1077 ! 1078 ALLOCATE(lon_rad(iml_dyn)) 1079 ALLOCATE(lon_ini(iml_dyn)) 1080 1081 IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 1082 lon_ini(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0 1083 ELSE 1084 lon_ini(:) = lon_dyn(:,1) 1085 ENDIF 1086 1087 ALLOCATE(lat_rad(jml_dyn)) 1088 ALLOCATE(lat_ini(jml_dyn)) 1089 1090 ALLOCATE(lev_dyn(llm_dyn)) 1091 1092 IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 1093 lat_ini(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0 1094 ELSE 1095 lat_ini(:) = lat_dyn(1,:) 1096 ENDIF 1097 ! 1098 1099 CALL conf_dat3d ( varname,iml_dyn, jml_dyn, llm_dyn, lon_ini, 1100 . lat_ini, levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d , 1101 , interbar ) 1102 1103 ALLOCATE(var_tmp2d(iml-1, jml)) 1104 ALLOCATE(var_tmp3d(iml, jml, llm_dyn)) 1105 ALLOCATE(ax(llm_dyn)) 1106 ALLOCATE(ay(llm_dyn)) 1107 ALLOCATE(yder(llm_dyn)) 1108 ALLOCATE(lind(llm_dyn)) 1109 ! 1110 1111 DO il=1,llm_dyn 1112 ! 1113 IF( interbar ) THEN 1114 IF( il.EQ.1 ) THEN 1115 WRITE(6,*) '-------------------------------------------------', 1116 ,'--------------' 1117 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 1118 , ' pour ', varname 1119 WRITE(6,*) '-------------------------------------------------', 1120 ,'--------------' 1121 ENDIF 1122 CALL inter_barxy ( iml_dyn, jml_dyn -1,lon_rad, lat_rad, 1123 , var_ana3d(:,:,il),iml-1, jml2, lon_in2, lat_in2,jml,var_tmp2d ) 1124 ELSE 1125 CALL grille_m(iml_dyn, jml_dyn, lon_rad, lat_rad, 1126 . var_ana3d(:,:,il), iml-1, jml, lon_in, lat_in, var_tmp2d ) 1127 ENDIF 1128 ! 1129 CALL gr_int_dyn(var_tmp2d, var_tmp3d(:,:,il), iml-1, jml) 1130 ! 1131 ENDDO 1132 ! 1133 DO il=1,llm_dyn 1134 lind(il) = llm_dyn-il+1 1135 ENDDO 1136 ! 1137 c 1138 c ... Pour l'interpolation verticale ,on interpole du haut de l'atmosphere 1139 c vers le sol ... 1140 c 1141 DO ij=1,jml 1142 DO ii=1,iml-1 1143 ! 1144 ax(:) = lev_dyn(lind(:)) 1145 ay(:) = var_tmp3d(ii, ij, lind(:)) 1146 ! 1147 1148 CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder) 1149 ! 1150 DO il=1,lml 1151 bx = pls_in(ii, ij, il) 1152 CALL SPLINT(ax, ay, yder, llm_dyn, bx, by) 1153 var3d(ii, ij, il) = by 1154 ENDDO 1155 ! 1156 ENDDO 1157 var3d(iml, ij, :) = var3d(1, ij, :) 1158 ENDDO 1159 1160 do il=1,lml 1161 call minmax(iml*jml,var3d(1,1,il),chmin,chmax) 1162 SELECTCASE(varname) 1163 CASE('U') 1164 WRITE(*,*) ' U min max l ',il,chmin,chmax 1165 CASE('V') 1166 WRITE(*,*) ' V min max l ',il,chmin,chmax 1167 CASE('TEMP') 1168 WRITE(*,*) ' TEMP min max l ',il,chmin,chmax 1169 CASE('R') 1170 WRITE(*,*) ' R min max l ',il,chmin,chmax 1171 END SELECT 1172 enddo 1173 1174 DEALLOCATE(lon_rad) 1175 DEALLOCATE(lon_ini) 1176 DEALLOCATE(lat_rad) 1177 DEALLOCATE(lat_ini) 1178 DEALLOCATE(lev_dyn) 1179 DEALLOCATE(var_tmp2d) 1180 DEALLOCATE(var_tmp3d) 1181 DEALLOCATE(ax) 1182 DEALLOCATE(ay) 1183 DEALLOCATE(yder) 1184 DEALLOCATE(lind) 1185 1186 ! 1187 RETURN 1188 ! 1189 END SUBROUTINE start_inter_3d 1190 ! 530 CHARACTER(LEN=25) :: title 531 CHARACTER(LEN=120) :: physfname 532 LOGICAL :: check=.TRUE. 533 REAL :: date, dt 534 INTEGER, DIMENSION(1) :: itau 535 INTEGER :: i, j 536 REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana, z 537 REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini 538 REAL, DIMENSION(:), ALLOCATABLE :: xppn, xpps 539 !------------------------------------------------------------------------------- 540 541 !--- KINETIC ENERGY 542 physfname = 'ECDYN.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0 543 IF(check) WRITE(lunout,*) 'Opening the surface analysis' 544 CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn) 545 IF(check) WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn 546 547 ALLOCATE(lat_dyn(iml_dyn,jml_dyn)) 548 ALLOCATE(lon_dyn(iml_dyn,jml_dyn)) 549 ALLOCATE(levdyn_ini(llm_dyn)) 550 CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, lon_dyn,lat_dyn,& 551 levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn) 552 553 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 554 ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn)) 555 lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad 556 lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad 557 558 ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn)) 559 560 !--- SURFACE GEOPOTENTIAL 561 title='Z' 562 ALLOCATE(z(iml,jml)) 563 CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana) 564 CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, & 565 var_ana, ibar) 566 CALL interp_startvar(title, ibar, .TRUE., & 567 iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & 568 lon_in, lat_in, lon_in2, lat_in2, z) 569 570 !--- SURFACE PRESSURE 571 title='SP' 572 ALLOCATE(psol_dyn(iml,jml)) 573 CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana) 574 CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, & 575 var_ana, ibar) 576 CALL interp_startvar(title, ibar, .TRUE., & 577 iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & 578 lon_in, lat_in, lon_in2, lat_in2, psol_dyn) 579 580 DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) 581 582 !--- ALLOCATION OF VARIABLES CREATED IN OR COMING FROM RESTART FILE 583 IF(.NOT.ALLOCATED(tsol)) THEN 584 CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 585 ELSE 586 IF(SIZE(tsol)/=SIZE(psol_dyn)) THEN 587 WRITE(lunout,*)'start_init_dyn :' 588 WRITE(lunout,*)'The temperature field we have does not have the right size' 589 STOP 590 END IF 591 END IF 592 593 IF(.NOT.ALLOCATED(phis)) THEN 594 CALL start_init_orog(iml,jml,lon_in,lat_in) 595 ELSE 596 IF(SIZE(phis)/=SIZE(psol_dyn)) THEN 597 WRITE(lunout,*)'start_init_dyn :' 598 WRITE(lunout,*)'The orography field we have does not have the right size' 599 STOP 600 END IF 601 END IF 602 603 !--- PSOL IS COMPUTED IN PASCALS 604 DO j = 1, jml 605 DO i = 1, iml-1 606 psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))/287.0/tsol(i,j)) 607 END DO 608 psol_dyn(iml,j) = psol_dyn(1,j) 609 END DO 610 DEALLOCATE(z) 611 612 ALLOCATE(xppn(iml-1),xpps(iml-1)) 613 DO i = 1, iml-1 614 xppn(i) = aire( i,1) * psol_dyn( i,1) 615 xpps(i) = aire( i,jml) * psol_dyn( i,jml) 616 END DO 617 DO i = 1, iml 618 psol_dyn(i,1 ) = SUM(xppn)/apoln 619 psol_dyn(i,jml) = SUM(xpps)/apols 620 END DO 621 DEALLOCATE(xppn,xpps) 622 623 RETURN 624 625 END SUBROUTINE start_init_dyn 626 ! 627 !------------------------------------------------------------------------------- 628 629 630 !------------------------------------------------------------------------------- 631 ! 632 SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,& 633 lat_in2, pls_in, var3d, ibar) 634 ! 635 !------------------------------------------------------------------------------- 636 ! Arguments: 637 CHARACTER(LEN=*), INTENT(IN) :: varname 638 INTEGER, INTENT(IN) :: iml, jml, lml 639 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 640 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 641 INTEGER, INTENT(IN) :: jml2 642 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 643 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 644 REAL, DIMENSION(iml,jml,lml), INTENT(IN) :: pls_in 645 REAL, DIMENSION(iml,jml,lml), INTENT(OUT) :: var3d 646 LOGICAL, INTENT(IN) :: ibar 647 !------------------------------------------------------------------------------- 648 ! Local variables: 649 #include "iniprint.h" 650 LOGICAL :: check=.TRUE. 651 REAL :: bx, by, chmin, chmax 652 INTEGER :: ii, ij, il 653 REAL, DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d 654 REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini 655 REAL, DIMENSION(:), ALLOCATABLE :: lev_dyn, ax, ay, yder 656 INTEGER, DIMENSION(:), ALLOCATABLE :: lind 657 !------------------------------------------------------------------------------- 658 IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D field.' 659 IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn 660 IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn 661 662 IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn)) 663 CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d) 664 665 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 666 ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn)) 667 lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad 668 lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad 669 670 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS 671 ALLOCATE(lon_rad(iml_dyn),lat_rad(jml_dyn),lev_dyn(llm_dyn)) 672 CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini, & 673 levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar) 674 DEALLOCATE(lon_ini,lat_ini) 675 676 !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro 677 ALLOCATE(var_tmp3d(iml,jml,llm_dyn)) 678 DO il=1,llm_dyn 679 CALL interp_startvar(varname, ibar, il==1, & 680 iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2, & 681 lon_in, lat_in, lon_in2, lat_in2, var_tmp3d(:,:,il)) 682 END DO 683 DEALLOCATE(lon_rad,lat_rad) 684 685 ALLOCATE(lind(llm_dyn)) 686 DO il=1,llm_dyn 687 lind(il) = llm_dyn-il+1 688 END DO 689 690 !--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND 691 ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn)) 692 DO ij=1,jml 693 DO ii=1,iml-1 694 ax(:)=lev_dyn(lind(:)) 695 ay(:)=var_tmp3d(ii,ij,lind(:)) 696 CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder) 697 DO il=1,lml 698 bx=pls_in(ii,ij,il) 699 CALL SPLINT(ax, ay, yder, llm_dyn, bx, by) 700 var3d(ii,ij,il) = by 701 END DO 702 END DO 703 var3d(iml,ij,:) = var3d(1,ij,:) 704 END DO 705 DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind) 706 707 DO il=1,lml 708 CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax) 709 WRITE(lunout,*)' '//TRIM(varname)//' min max l ',il,chmin,chmax 710 END DO 711 712 RETURN 713 714 END SUBROUTINE start_inter_3d 715 ! 716 !------------------------------------------------------------------------------- 717 718 719 !------------------------------------------------------------------------------- 720 ! 721 SUBROUTINE interp_startvar(vname, ibar, ibeg, ii, jj, lon, lat, vari, & 722 i1, j1, j2, lon1, lat1, lon2, lat2, varo) 723 ! 724 !------------------------------------------------------------------------------- 725 ! Arguments: 726 CHARACTER(LEN=*), INTENT(IN) :: vname 727 LOGICAL, INTENT(IN) :: ibar, ibeg 728 INTEGER, INTENT(IN) :: ii, jj 729 REAL, DIMENSION(ii), INTENT(IN) :: lon 730 REAL, DIMENSION(jj), INTENT(IN) :: lat 731 REAL, DIMENSION(ii,jj), INTENT(IN) :: vari 732 INTEGER, INTENT(IN) :: i1, j1, j2 733 REAL, DIMENSION(i1), INTENT(IN) :: lon1 734 REAL, DIMENSION(j1), INTENT(IN) :: lat1 735 REAL, DIMENSION(i1), INTENT(IN) :: lon2 736 REAL, DIMENSION(j2), INTENT(IN) :: lat2 737 REAL, DIMENSION(i1,j1), INTENT(OUT) :: varo 738 !------------------------------------------------------------------------------- 739 ! Local variables: 740 #include "iniprint.h" 741 REAL, DIMENSION(i1-1,j1) :: vtmp 742 !------------------------------------------------------------------------------- 743 IF(ibar) THEN 744 IF(ibeg) THEN 745 WRITE(lunout,*) & 746 '---------------------------------------------------------------' 747 WRITE(lunout,*) & 748 '$$$ Utilisation de l interpolation barycentrique pour '//TRIM(vname)//' $$$' 749 WRITE(lunout,*) & 750 '---------------------------------------------------------------' 751 END IF 752 CALL inter_barxy(ii, jj-1, lon, lat, vari, i1-1, j2, lon2, lat2, j1, vtmp) 753 ELSE 754 CALL grille_m (ii, jj, lon, lat, vari, i1-1, j1, lon1, lat1, vtmp) 755 END IF 756 CALL gr_int_dyn(vtmp, varo, i1-1, j1) 757 758 END SUBROUTINE interp_startvar 759 ! 760 !------------------------------------------------------------------------------- 761 1191 762 #endif 1192 763 ! of #ifdef CPP_EARTH 1193 END MODULE startvar 764 765 END MODULE startvar 766 ! 767 !*******************************************************************************
Note: See TracChangeset
for help on using the changeset viewer.
