Changeset 2594 for LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd
- Timestamp:
- Jul 18, 2016, 9:41:10 PM (8 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 3 edited
-
. (modified) (1 prop)
-
libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90 (modified) (4 diffs)
-
libf/dynphy_lonlat/phylmd/limit_netcdf.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 2547-2567,2569,2571-2574,2576-2589
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90
r2471 r2594 12 12 prad,pg,pr,pcpp,iflag_phys) 13 13 USE dimphy, ONLY: init_dimphy 14 USE mod_grid_phy_lmdz, ONLY: klon_glo, & ! number of atmospheric columns (on full grid) 15 regular_lonlat, & ! regular longitude-latitude grid type 16 nbp_lon, nbp_lat, nbp_lev 17 USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid) 18 klon_omp_begin, & ! start index of local omp subgrid 19 klon_omp_end, & ! end index of local omp subgrid 20 klon_mpi_begin ! start indes of columns (on local mpi grid) 21 USE geometry_mod, ONLY : init_geometry 14 USE inigeomphy_mod, ONLY: inigeomphy 15 USE mod_grid_phy_lmdz, ONLY: klon_glo ! number of atmospheric columns (on full grid) 16 USE mod_phys_lmdz_para, ONLY: klon_omp ! number of columns (on local omp grid) 22 17 USE vertical_layers_mod, ONLY : init_vertical_layers 23 18 USE infotrac, ONLY: nqtot,nqo,nbtr,tname,ttext,type_trac,& … … 39 34 USE phystokenc_mod, ONLY: init_phystokenc 40 35 USE phyaqua_mod, ONLY: iniaqua 41 USE physics_distribution_mod, ONLY : init_physics_distribution42 USE regular_lonlat_mod, ONLY : init_regular_lonlat, &43 east, west, north, south, &44 north_east, north_west, &45 south_west, south_east46 USE mod_interface_dyn_phys, ONLY : init_interface_dyn_phys47 36 #ifdef INCA 48 37 USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic … … 92 81 CHARACTER (LEN=20) :: modname = 'iniphysiq' 93 82 CHARACTER (LEN=80) :: abort_message 94 REAL :: total_area_phy, total_area_dyn95 83 96 ! boundaries, on global grid97 REAL,ALLOCATABLE :: boundslon_reg(:,:)98 REAL,ALLOCATABLE :: boundslat_reg(:,:)99 100 ! global array, on full physics grid:101 REAL,ALLOCATABLE :: latfi_glo(:)102 REAL,ALLOCATABLE :: lonfi_glo(:)103 REAL,ALLOCATABLE :: cufi_glo(:)104 REAL,ALLOCATABLE :: cvfi_glo(:)105 REAL,ALLOCATABLE :: airefi_glo(:)106 REAL,ALLOCATABLE :: boundslonfi_glo(:,:)107 REAL,ALLOCATABLE :: boundslatfi_glo(:,:)108 109 ! local arrays, on given MPI/OpenMP domain:110 REAL,ALLOCATABLE,SAVE :: latfi(:)111 REAL,ALLOCATABLE,SAVE :: lonfi(:)112 REAL,ALLOCATABLE,SAVE :: cufi(:)113 REAL,ALLOCATABLE,SAVE :: cvfi(:)114 REAL,ALLOCATABLE,SAVE :: airefi(:)115 REAL,ALLOCATABLE,SAVE :: boundslonfi(:,:)116 REAL,ALLOCATABLE,SAVE :: boundslatfi(:,:)117 !$OMP THREADPRIVATE (latfi,lonfi,cufi,cvfi,airefi,boundslonfi,boundslatfi)118 84 119 85 #ifndef CPP_PARA … … 123 89 #endif 124 90 125 ! Initialize Physics distibution and parameters and interface with dynamics126 IF (ii*jj>1) THEN ! general 3D case127 CALL init_physics_distribution(regular_lonlat,4, &128 nbp,ii,jj+1,nlayer,communicator)129 ELSE ! For 1D model130 CALL init_physics_distribution(regular_lonlat,4, &131 1,1,1,nlayer,communicator)132 ENDIF 133 CALL init_interface_dyn_phys91 ! --> initialize physics distribution, global fields and geometry 92 ! (i.e. things in phy_common or dynphy_lonlat) 93 CALL inigeomphy(ii,jj,nlayer, & 94 nbp, communicator, & 95 rlatu,rlatv, & 96 rlonu,rlonv, & 97 aire,cu,cv) 98 99 ! --> now initialize things specific to the phylmd physics package 134 100 135 ! init regular global longitude-latitude grid points and boundaries136 ALLOCATE(boundslon_reg(ii,2))137 ALLOCATE(boundslat_reg(jj+1,2))138 139 DO i=1,ii140 boundslon_reg(i,east)=rlonu(i)141 boundslon_reg(i,west)=rlonu(i+1)142 ENDDO143 144 boundslat_reg(1,north)= PI/2145 boundslat_reg(1,south)= rlatv(1)146 DO j=2,jj147 boundslat_reg(j,north)=rlatv(j-1)148 boundslat_reg(j,south)=rlatv(j)149 ENDDO150 boundslat_reg(jj+1,north)= rlatv(jj)151 boundslat_reg(jj+1,south)= -PI/2152 153 ! Write values in module regular_lonlat_mod154 CALL init_regular_lonlat(ii,jj+1, rlonv(1:ii), rlatu, &155 boundslon_reg, boundslat_reg)156 157 ! Generate global arrays on full physics grid158 ALLOCATE(latfi_glo(klon_glo),lonfi_glo(klon_glo))159 ALLOCATE(cufi_glo(klon_glo),cvfi_glo(klon_glo))160 ALLOCATE(airefi_glo(klon_glo))161 ALLOCATE(boundslonfi_glo(klon_glo,4))162 ALLOCATE(boundslatfi_glo(klon_glo,4))163 164 165 IF (klon_glo>1) THEN ! general case166 ! North pole167 latfi_glo(1)=rlatu(1)168 lonfi_glo(1)=0.169 cufi_glo(1) = cu(1)170 cvfi_glo(1) = cv(1)171 boundslonfi_glo(1,north_east)=0172 boundslatfi_glo(1,north_east)=PI/2173 boundslonfi_glo(1,north_west)=2*PI174 boundslatfi_glo(1,north_west)=PI/2175 boundslonfi_glo(1,south_west)=2*PI176 boundslatfi_glo(1,south_west)=rlatv(1)177 boundslonfi_glo(1,south_east)=0178 boundslatfi_glo(1,south_east)=rlatv(1)179 DO j=2,jj180 DO i=1,ii181 k=(j-2)*ii+1+i182 latfi_glo(k)= rlatu(j)183 lonfi_glo(k)= rlonv(i)184 cufi_glo(k) = cu((j-1)*(ii+1)+i)185 cvfi_glo(k) = cv((j-1)*(ii+1)+i)186 boundslonfi_glo(k,north_east)=rlonu(i)187 boundslatfi_glo(k,north_east)=rlatv(j-1)188 boundslonfi_glo(k,north_west)=rlonu(i+1)189 boundslatfi_glo(k,north_west)=rlatv(j-1)190 boundslonfi_glo(k,south_west)=rlonu(i+1)191 boundslatfi_glo(k,south_west)=rlatv(j)192 boundslonfi_glo(k,south_east)=rlonu(i)193 boundslatfi_glo(k,south_east)=rlatv(j)194 ENDDO195 ENDDO196 ! South pole197 latfi_glo(klon_glo)= rlatu(jj+1)198 lonfi_glo(klon_glo)= 0.199 cufi_glo(klon_glo) = cu((ii+1)*jj+1)200 cvfi_glo(klon_glo) = cv((ii+1)*jj-ii)201 boundslonfi_glo(klon_glo,north_east)= 0202 boundslatfi_glo(klon_glo,north_east)= rlatv(jj)203 boundslonfi_glo(klon_glo,north_west)= 2*PI204 boundslatfi_glo(klon_glo,north_west)= rlatv(jj)205 boundslonfi_glo(klon_glo,south_west)= 2*PI206 boundslatfi_glo(klon_glo,south_west)= -PI/2207 boundslonfi_glo(klon_glo,south_east)= 0208 boundslatfi_glo(klon_glo,south_east)= -Pi/2209 210 ! build airefi_glo(), mesh area on physics grid211 CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi_glo)212 ! Poles are single points on physics grid213 airefi_glo(1)=sum(aire(1:ii,1))214 airefi_glo(klon_glo)=sum(aire(1:ii,jj+1))215 216 ! Sanity check: do total planet area match between physics and dynamics?217 total_area_dyn=sum(aire(1:ii,1:jj+1))218 total_area_phy=sum(airefi_glo(1:klon_glo))219 IF (total_area_dyn/=total_area_phy) THEN220 WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'221 WRITE (lunout, *) ' in the dynamics total_area_dyn=', total_area_dyn222 WRITE (lunout, *) ' but in the physics total_area_phy=', total_area_phy223 IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN224 ! stop here if the relative difference is more than 0.001%225 abort_message = 'planet total surface discrepancy'226 CALL abort_gcm(modname, abort_message, 1)227 ENDIF228 ENDIF229 ELSE ! klon_glo==1, running the 1D model230 ! just copy over input values231 latfi_glo(1)=rlatu(1)232 lonfi_glo(1)=rlonv(1)233 cufi_glo(1)=cu(1)234 cvfi_glo(1)=cv(1)235 airefi_glo(1)=aire(1,1)236 boundslonfi_glo(1,north_east)=rlonu(1)237 boundslatfi_glo(1,north_east)=PI/2238 boundslonfi_glo(1,north_west)=rlonu(2)239 boundslatfi_glo(1,north_west)=PI/2240 boundslonfi_glo(1,south_west)=rlonu(2)241 boundslatfi_glo(1,south_west)=rlatv(1)242 boundslonfi_glo(1,south_east)=rlonu(1)243 boundslatfi_glo(1,south_east)=rlatv(1)244 ENDIF ! of IF (klon_glo>1)245 246 101 !$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/) 247 ! Now generate local lon/lat/cu/cv/area/bounds arrays248 ALLOCATE(latfi(klon_omp),lonfi(klon_omp),cufi(klon_omp),cvfi(klon_omp))249 ALLOCATE(airefi(klon_omp))250 ALLOCATE(boundslonfi(klon_omp,4))251 ALLOCATE(boundslatfi(klon_omp,4))252 253 254 offset = klon_mpi_begin - 1255 airefi(1:klon_omp) = airefi_glo(offset+klon_omp_begin:offset+klon_omp_end)256 cufi(1:klon_omp) = cufi_glo(offset+klon_omp_begin:offset+klon_omp_end)257 cvfi(1:klon_omp) = cvfi_glo(offset+klon_omp_begin:offset+klon_omp_end)258 lonfi(1:klon_omp) = lonfi_glo(offset+klon_omp_begin:offset+klon_omp_end)259 latfi(1:klon_omp) = latfi_glo(offset+klon_omp_begin:offset+klon_omp_end)260 boundslonfi(1:klon_omp,:) = boundslonfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)261 boundslatfi(1:klon_omp,:) = boundslatfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)262 263 ! copy over local grid longitudes and latitudes264 CALL init_geometry(klon_omp,lonfi,latfi,boundslonfi,boundslatfi, &265 airefi,cufi,cvfi)266 102 267 103 ! copy over preff , ap(), bp(), etc -
LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/limit_netcdf.F90
r2542 r2594 117 117 ALLOCATE(pctsrf(klon,nbsrf)) 118 118 CALL start_init_subsurf(.FALSE.) 119 !--- TO MATCH EXACTLY WHAT WOULD BE DONE IN etat0phys_netcdf 120 WHERE( masque(:,:)<EPSFRA) masque(:,:)=0. 121 WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1. 119 122 END IF 120 123 … … 336 339 REAL, ALLOCATABLE :: champan(:,:,:) 337 340 !--- input files 341 CHARACTER(LEN=20) :: fnam_m, fnam_p ! previous/next files names 338 342 CHARACTER(LEN=20) :: cal_in ! calendar 339 343 CHARACTER(LEN=20) :: unit_sic ! attribute unit in sea-ice file … … 345 349 LOGICAL :: extrp ! flag for extrapolation 346 350 REAL :: chmin, chmax 347 INTEGER ierr 351 INTEGER ierr, idx 348 352 integer n_extrap ! number of extrapolated points 349 353 logical skip … … 360 364 END SELECT 361 365 extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag 366 idx=INDEX(fnam,'.nc')-1 362 367 363 368 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ----------------------------- … … 393 398 !--- Time (variable is not needed - it is rebuilt - but calendar is) 394 399 CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam) 395 ALLOCATE(timeyear( MAX(lmdep,14)))400 ALLOCATE(timeyear(lmdep+2)) 396 401 CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam) 397 402 cal_in=' ' … … 412 417 IF(lmdep==12) THEN 413 418 timeyear=mid_month(anneeref, cal_in, ndays_in) 414 CALL msg(0,'Monthly periodic input file (perpetual run).') 415 ELSE IF(lmdep==14) THEN 416 timeyear=mid_month(anneeref, cal_in, ndays_in) 417 CALL msg(0,'Monthly 14-records input file (interannual run).') 419 CALL msg(0,'Monthly input file(s) for '//TRIM(title)//'.') 418 420 ELSE IF(lmdep==ndays_in) THEN 419 timeyear=[(REAL(k)-0.5,k= 1,ndays_in)]421 timeyear=[(REAL(k)-0.5,k=0,ndays_in+1)] 420 422 CALL msg(0,'Daily input file (no time interpolation).') 421 423 ELSE 422 424 WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep, & 423 ' records, 12/ 14/',ndays_in,' (periodic/interannual/daily) needed'425 ' records, 12/',ndays_in,' (monthly/daily needed).' 424 426 CALL abort_physic('mid_months',TRIM(mess),1) 425 427 END IF 426 428 427 429 !--- GETTING THE FIELD AND INTERPOLATING IT ---------------------------------- 428 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, MAX(lmdep,14)))430 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep+2)) 429 431 IF(extrp) ALLOCATE(work(imdep, jmdep)) 430 432 CALL msg(5,'') … … 446 448 WHERE(NINT(mask)/=1) champint=0.001 447 449 END IF 448 !--- Special case for periodic input file: index shifted 449 ll = l; IF(lmdep==12) ll = l + 1 450 champtime(:, :, ll)=champint 450 champtime(:, :, l+1)=champint 451 451 END DO 452 IF(lmdep==12) THEN453 champtime(:,:, 1)=champtime(:,:,13)454 champtime(:,:,14)=champtime(:,:, 2)455 END IF456 452 CALL ncerr(NF90_CLOSE(ncid), fnam) 457 453 454 !--- FIRST RECORD: LAST ONE OF PREVIOUS YEAR (CURRENT YEAR IF UNAVAILABLE) 455 fnam_m=fnam(1:idx)//'_m.nc' 456 IF(NF90_OPEN(fnam_m,NF90_NOWRITE,ncid)==NF90_NOERR) THEN 457 CALL msg(0,'Reading previous year file ("'//TRIM(fnam_m)//'") last record for '//TRIM(title)) 458 CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_m) 459 CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam_m) 460 CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), len=l), fnam_m) 461 CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam_m) 462 CALL ncerr(NF90_CLOSE(ncid), fnam_m) 463 CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.) 464 IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work) 465 IF(mode=='RUG') champ=LOG(champ) 466 CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint) 467 IF(mode=='RUG') THEN 468 champint=EXP(champint) 469 WHERE(NINT(mask)/=1) champint=0.001 470 END IF 471 champtime(:, :, 1)=champint 472 ELSE 473 CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") last record for '//TRIM(title)) 474 champtime(:, :, 1)=champtime(:, :, lmdep+1) 475 END IF 476 477 !--- LAST RECORD: FIRST ONE OF NEXT YEAR (CURRENT YEAR IF UNAVAILABLE) 478 fnam_p=fnam(1:idx)//'_p.nc' 479 IF(NF90_OPEN(fnam_p,NF90_NOWRITE,ncid)==NF90_NOERR) THEN 480 CALL msg(0,'Reading previous year file ("'//TRIM(fnam_p)//'") first record for '//TRIM(title)) 481 CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_p) 482 CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,1],[imdep,jmdep,1]),fnam_p) 483 CALL ncerr(NF90_CLOSE(ncid), fnam_p) 484 CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.) 485 IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work) 486 IF(mode=='RUG') champ=LOG(champ) 487 CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint) 488 IF(mode=='RUG') THEN 489 champint=EXP(champint) 490 WHERE(NINT(mask)/=1) champint=0.001 491 END IF 492 champtime(:, :, lmdep+2)=champint 493 ELSE 494 CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") first record for '//TRIM(title)) 495 champtime(:, :, lmdep+2)=champtime(:, :, 2) 496 END IF 458 497 DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ) 459 498 IF(extrp) DEALLOCATE(work) … … 477 516 END IF 478 517 END IF 479 ALLOCATE(yder( MAX(lmdep,14)), champan(iip1, jjp1, ndays))518 ALLOCATE(yder(lmdep+2), champan(iip1, jjp1, ndays)) 480 519 IF(lmdep==ndays_in) THEN 481 520 champan(1:iim,:,:)=champtime … … 546 585 ! 547 586 !------------------------------------------------------------------------------- 587 USE grid_noro_m, ONLY: grid_noro0 548 588 IMPLICIT NONE 549 589 !=============================================================================== … … 599 639 600 640 END SUBROUTINE start_init_orog0 601 !602 !-------------------------------------------------------------------------------603 604 605 !-------------------------------------------------------------------------------606 !607 SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)608 !609 !===============================================================================610 ! Purpose: Extracted from grid_noro to provide geopotential height for dynamics611 ! without any call to physics subroutines.612 !===============================================================================613 IMPLICIT NONE614 !-------------------------------------------------------------------------------615 ! Arguments:616 REAL, INTENT(IN) :: xd(:), yd(:) !--- INPUT COORDINATES (imdp) (jmdp)617 REAL, INTENT(IN) :: zd(:,:) !--- INPUT FIELD (imdp,jmdp)618 REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar)619 REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar)620 REAL, INTENT(INOUT):: mask(:,:) !--- MASK (imar+1,jmar)621 !-------------------------------------------------------------------------------622 ! Local variables:623 CHARACTER(LEN=256) :: modname="grid_noro0"624 REAL, ALLOCATABLE :: xusn(:), yusn(:) ! dim (imdp+2*iext) (jmdp+2)625 REAL, ALLOCATABLE :: zusn(:,:) ! dim (imdp+2*iext,jmdp+2)626 REAL, ALLOCATABLE :: weight(:,:) ! dim (imar+1,jmar)627 REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:) ! dim (imar+1,jmar)628 REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)629 REAL, ALLOCATABLE :: a(:), b(:) ! dim (imax)630 REAL, ALLOCATABLE :: c(:), d(:) ! dim (jmax)631 LOGICAL :: masque_lu632 INTEGER :: i, ii, imdp, imar, iext633 INTEGER :: j, jj, jmdp, jmar, nn634 REAL :: xpi, zlenx, weighx, xincr, zbordnor, zmeanor, zweinor, zbordest635 REAL :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue636 !-------------------------------------------------------------------------------637 imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")638 jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")639 imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1640 jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")641 IF(imar/=iim) CALL abort_gcm(TRIM(modname),'imar/=iim' ,1)642 IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)643 iext=imdp/10644 xpi = ACOS(-1.)645 rad = 6371229.646 647 !--- ARE WE USING A READ MASK ?648 masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0649 WRITE(lunout,*)'Masque lu: ',masque_lu650 651 !--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:652 ALLOCATE(xusn(imdp+2*iext))653 xusn(1 +iext:imdp +iext)=xd(:)654 xusn(1 : iext)=xd(1+imdp-iext:imdp)-2.*xpi655 xusn(1+imdp+iext:imdp+2*iext)=xd(1 :iext)+2.*xpi656 657 ALLOCATE(yusn(jmdp+2))658 yusn(1 )=yd(1) +(yd(1) -yd(2))659 yusn(2:jmdp+1)=yd(:)660 yusn( jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))661 662 ALLOCATE(zusn(imdp+2*iext,jmdp+2))663 zusn(1 +iext:imdp +iext,2:jmdp+1)=zd (: , :)664 zusn(1 : iext,2:jmdp+1)=zd (imdp-iext+1:imdp , :)665 zusn(1+imdp +iext:imdp+2*iext,2:jmdp+1)=zd (1:iext , :)666 zusn(1 :imdp/2+iext, 1)=zusn(1+imdp/2:imdp +iext, 2)667 zusn(1+imdp/2+iext:imdp+2*iext, 1)=zusn(1 :imdp/2+iext, 2)668 zusn(1 :imdp/2+iext, jmdp+2)=zusn(1+imdp/2:imdp +iext,jmdp+1)669 zusn(1+imdp/2+iext:imdp+2*iext, jmdp+2)=zusn(1 :imdp/2+iext,jmdp+1)670 671 !--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)672 ALLOCATE(a(imar+1),b(imar+1))673 b(1:imar)=(x(1:imar )+ x(2:imar+1))/2.0674 b(imar+1)= x( imar+1)+(x( imar+1)-x(imar))/2.0675 a(1)=x(1)-(x(2)-x(1))/2.0676 a(2:imar+1)= b(1:imar)677 678 ALLOCATE(c(jmar),d(jmar))679 d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0680 d( jmar )= y( jmar )+(y( jmar)-y(jmar-1))/2.0681 c(1)=y(1)-(y(2)-y(1))/2.0682 c(2:jmar)=d(1:jmar-1)683 684 !--- INITIALIZATIONS:685 ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0686 ALLOCATE(zmea (imar+1,jmar)); zmea (:,:)= 0.0687 688 !--- SUMMATION OVER GRIDPOINT AREA689 zleny=xpi/REAL(jmdp)*rad690 xincr=xpi/REAL(jmdp)/2.691 ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.692 ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.693 DO ii = 1, imar+1694 DO jj = 1, jmar695 DO j = 2,jmdp+1696 zlenx =zleny *COS(yusn(j))697 zbordnor=(xincr+c(jj)-yusn(j))*rad698 zbordsud=(xincr-d(jj)+yusn(j))*rad699 weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))700 IF(weighy/=0) THEN701 DO i = 2, imdp+2*iext-1702 zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))703 zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))704 weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))705 IF(weighx/=0)THEN706 num_tot(ii,jj)=num_tot(ii,jj)+1.0707 IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0708 weight(ii,jj)=weight(ii,jj)+weighx*weighy709 zmea (ii,jj)=zmea (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN710 END IF711 END DO712 END IF713 END DO714 END DO715 END DO716 717 !--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME718 IF(.NOT.masque_lu) THEN719 WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)720 END IF721 nn=COUNT(weight(:,1:jmar-1)==0.0)722 IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn723 WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)724 725 !--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)726 ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0727 WHERE(mask>=0.1) mask_tmp = 1.728 WHERE(weight(:,:)/=0.0)729 zphi(:,:)=mask_tmp(:,:)*zmea(:,:)730 zmea(:,:)=mask_tmp(:,:)*zmea(:,:)731 END WHERE732 WRITE(lunout,*)' MEAN ORO:' ,MAXVAL(zmea)733 734 !--- Values at poles735 zphi(imar+1,:)=zphi(1,:)736 737 zweinor=SUM(weight(1:imar, 1),DIM=1)738 zweisud=SUM(weight(1:imar,jmar),DIM=1)739 zmeanor=SUM(weight(1:imar, 1)*zmea(1:imar, 1),DIM=1)740 zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)741 zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud742 743 END SUBROUTINE grid_noro0744 641 ! 745 642 !-------------------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
