Changeset 7796 for branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/globgrd.f90
- Timestamp:
- 2022-11-08T16:24:46+01:00 (20 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/globgrd.f90
r7261 r7796 80 80 ! 81 81 INTEGER(i_std) :: iret, ndims, nvars, nb_atts, id_unlim, iv, lll 82 INTEGER(i_std) :: iindex_init, jindex_init, iindex_end, jindex_end 82 83 INTEGER(i_std) :: iim_full, jjm_full, nbland_full 83 84 CHARACTER(LEN=20) :: axname, varname 84 85 CHARACTER(LEN=120) :: tmpfile 85 86 REAL(r_std), DIMENSION(2) :: loczoom_lon, loczoom_lat 86 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask 87 INTEGER(i_std), DIMENSION(2) :: tmp_lon_ind1, tmp_lat_ind1, tmp_lon_ind2, tmp_lat_ind2 88 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask,zoom_mask, lat, lon 87 89 ! 88 90 ! Set default values against which we can test … … 92 94 ! 93 95 ! Verify the grid file name 96 ! if forcing_file exists and if grid file exists: tmpfile = grid file 97 ! if forcing_file exits and grid file not exists: tmpfile = forcing_file 98 ! if forcing_file not exists: tmp_file= grid file 94 99 ! 95 IF ( INDEX(filename,"NONE") >= 1) THEN100 IF ( PRESENT(forcingfile) ) THEN 96 101 is_forcing_file=.TRUE. 97 IF ( PRESENT(forcingfile)) THEN102 IF ( INDEX(filename,"NONE") >= 1 ) THEN 98 103 tmpfile=forcingfile(1) 99 104 ELSE 100 CALL ipslerr (3,'globgrd_getdomsz',"Error No grid file provided :",tmpfile, " ")105 tmpfile=filename 101 106 ENDIF 102 107 ELSE … … 105 110 ENDIF 106 111 ! 107 ! Verify that the zo mmis provided. Else choose the entire globe112 ! Verify that the zoomed region is provided. Else choose the entire globe 108 113 ! 109 114 IF ( PRESENT(zoom_lon) .AND. PRESENT(zoom_lat) ) THEN … … 144 149 !! Coordinate variables used by WRF. 145 150 CASE("west_east") 146 iim = lll151 iim_full = lll 147 152 model_guess = "WRF" 148 153 CASE("south_north") 149 jjm = lll154 jjm_full = lll 150 155 model_guess = "WRF" 151 156 ! 152 157 !! Variables used in WFDEI 153 158 CASE("lon") 154 iim = lll159 iim_full = lll 155 160 model_guess = "regular" 156 161 CASE("lat") 157 jjm = lll162 jjm_full = lll 158 163 model_guess = "regular" 159 164 CASE("nbland") 160 nbland = lll165 nbland_full = lll 161 166 ! 162 167 !! Variables used by CRU-NCEP 163 168 CASE("nav_lon") 164 iim = lll169 iim_full = lll 165 170 model_guess = "regular" 166 171 CASE("nav_lat") 167 jjm = lll172 jjm_full = lll 168 173 model_guess = "regular" 169 174 CASE("land") 170 nbland = lll 175 nbland_full = lll 176 177 178 !! Variables used by CRU-JRA (v2.2.2) 179 CASE("longitude") 180 iim_full = lll 181 model_guess = "regular" 182 CASE("latitude") 183 jjm_full = lll 184 model_guess = "regular" 171 185 END SELECT 172 186 ENDDO 173 187 ! 174 ! If we have a WRF file we need to count the number of land points .188 ! If we have a WRF file we need to count the number of land points, define iim jjm and mask 175 189 ! 176 190 IF ( model_guess == "WRF" ) THEN 177 191 178 ALLOCATE(mask(iim,jjm))192 IF ( .NOT. ALLOCATED(mask) ) ALLOCATE(mask(iim_full,jjm_full)) 179 193 180 194 varname = "LANDMASK" … … 186 200 ENDIF 187 201 188 nbland = COUNT(mask > 0.5) 189 190 ENDIF 191 ! 192 ! 193 ! If we are in the case of a forcing file a few functions from forcing_tools need to be called 202 nbland_full = COUNT(mask > 0.5) 203 204 205 IF ( .NOT. ALLOCATED(lat)) ALLOCATE(lat(iim_full,jjm_full)) 206 IF ( .NOT. ALLOCATED(lon)) ALLOCATE(lon(iim_full,jjm_full)) 207 208 varname = "XLONG_M" 209 iret = NF90_INQ_VARID (fid, varname, iv) 210 IF (iret /= NF90_NOERR) THEN 211 CALL ipslerr (3,'globgrd_getdomsz',"Could not find variable ", varname," ") 212 ELSE 213 iret = NF90_GET_VAR (fid,iv,lon) 214 ENDIF 215 216 varname = "XLAT_M" 217 iret = NF90_INQ_VARID (fid, varname, iv) 218 IF (iret /= NF90_NOERR) THEN 219 CALL ipslerr (3,'globgrd_getdomsz',"Could not find variable ", varname," ") 220 ELSE 221 iret = NF90_GET_VAR (fid,iv,lat) 222 ENDIF 223 224 !define nbland for full-region simulation in case of need 225 nbland = nbland_full 226 227 ENDIF 228 ! 229 ! 230 ! If we are in the case of a forcing file and model_guess is regular grid, 231 ! then a few functions from forcing_tools need to be called 194 232 ! so that the file is analysed with the tools of the forcing module. 195 233 ! 196 IF ( is_forcing_file ) THEN 234 ! predefine nbland, iim, jjm in the case is_forcing_file=false, or full grid simulation for wrf 235 iim = iim_full 236 jjm = jjm_full 237 IF ( is_forcing_file .AND. model_guess .EQ. "regular") THEN 197 238 ! 198 239 ! Because we are re-using routines from the forcing module, we have to … … 204 245 ENDIF 205 246 ! 206 ! This has to be a regular grid. A more clever clasification of files will be needed.207 model_guess = "regular"208 209 247 ! Set last argument closefile=.FALSE. as the forcing file has been closed here above. 210 248 ! This will also induce that dump_mask=.FALSE. in forcing_getglogrid and the … … 213 251 WRITE(*,*) forcingfile, "Forcing file with dimensions : ", iim_full, jjm_full, nbland_full 214 252 ! 215 CALL forcing_zoomgrid (loczoom_lon, loczoom_lat, forcingfile(1), .TRUE.)253 CALL forcing_zoomgrid (loczoom_lon, loczoom_lat, forcingfile(1), model_guess, .TRUE.) 216 254 ! 217 255 CALL forcing_givegridsize (iim, jjm, nbland) 218 256 ! 257 258 ELSE IF ( is_forcing_file .AND. model_guess .EQ. "WRF") THEN 259 ! 260 ! if forcing_file exists and model_guess is wrf, and zoomed domain is asked, 261 ! then we define the domain size according to the zoom 262 ! if the full domain is asked, then we do not define the domain size again here. 263 ! 264 ! get the beginning and endding index in the case of zoomed region, for WRF grids (XW) 265 ! Not to close the grid file here, to be used by globgrd_getgrid 266 ! 267 !IF ( PRESENT(zoom_lon) .AND. PRESENT(zoom_lat) ) THEN 268 IF ( PRESENT(zoom_lon) .OR. PRESENT(zoom_lat) ) THEN 269 270 ! to get new iim, jjm and nbland for the zoomed region 271 tmp_lon_ind1 = MINLOC(ABS(lon(:,:)-loczoom_lon(1))) 272 tmp_lat_ind1 = MINLOC(ABS(lat(:,:)-loczoom_lat(1))) 273 tmp_lon_ind2 = MINLOC(ABS(lon(:,:)-loczoom_lon(2))) 274 tmp_lat_ind2 = MINLOC(ABS(lat(:,:)-loczoom_lat(2))) 275 ! 276 ! the zoomed region 277 ! lon1, lat2------- lon2,lat2 278 ! | | 279 ! lon1, lat1------- lon2,lat1 280 ! 281 iindex_init = tmp_lon_ind1(1) 282 jindex_init= tmp_lat_ind1(2) 283 284 iindex_end= tmp_lon_ind2(1) 285 jindex_end= tmp_lat_ind2(2) 286 287 iim = iindex_end - iindex_init + 1 288 jjm = jindex_end - jindex_init + 1 289 290 291 IF ( .NOT. ALLOCATED(zoom_mask) ) ALLOCATE(zoom_mask(iim,jjm)) 292 zoom_mask = mask(iindex_init:iindex_end, jindex_init:jindex_end) 293 nbland = COUNT(zoom_mask > 0.5) 294 IF ( ALLOCATED(zoom_mask) ) DEALLOCATE(zoom_mask) 295 296 ENDIF 297 298 IF ( ALLOCATED(lat) ) DEALLOCATE(lat) 299 IF ( ALLOCATED(lon) ) DEALLOCATE(lon) 300 IF ( ALLOCATED(mask) ) DEALLOCATE(mask) 301 219 302 ENDIF 220 303 ! … … 250 333 !--------------------------------------------------------------------- 251 334 SUBROUTINE globgrd_getgrid(fid, iim, jjm, nbland, model_guess, lon, lat, mask, area, corners, & 252 & lindex, contfrac, calendar )335 & lindex, contfrac, calendar, zoom_lon, zoom_lat) 253 336 ! 254 337 ! … … 262 345 INTEGER(i_std), INTENT(in) :: iim, jjm, nbland 263 346 CHARACTER(LEN=*), INTENT(in) :: model_guess 347 REAL(r_std), DIMENSION(2), INTENT(in), OPTIONAL :: zoom_lon, zoom_lat 264 348 ! 265 349 ! OUTPUT … … 275 359 CASE("WRF") 276 360 CALL globgrd_getwrf(fid, iim, jjm, nbland, lon, lat, mask, area, corners, & 277 & lindex, contfrac, calendar )361 & lindex, contfrac, calendar, zoom_lon, zoom_lat) 278 362 CASE("regular") 279 363 IF ( .NOT. is_forcing_file ) THEN … … 482 566 !! 483 567 SUBROUTINE globgrd_getwrf(fid, iim, jjm, nbland, lon, lat, mask, area, corners, & 484 & lindex, contfrac, calendar )568 & lindex, contfrac, calendar, zoom_lon, zoom_lat) 485 569 ! 486 570 USE defprec … … 491 575 INTEGER(i_std), INTENT(in) :: fid 492 576 INTEGER(i_std), INTENT(in) :: iim, jjm, nbland 577 REAL(r_std), DIMENSION(2), INTENT(in), OPTIONAL :: zoom_lon, zoom_lat 493 578 ! 494 579 ! OUTPUT … … 503 588 ! 504 589 INTEGER(i_std) :: i, ip, jp, k, iret, iv, nvars, varndim 505 INTEGER(i_std) :: imm1, imp1 590 INTEGER(i_std),DIMENSION(2) :: tmp_lon_ind1, tmp_lat_ind1, tmp_lon_ind2, tmp_lat_ind2 591 REAL(r_std), DIMENSION(2) :: loczoom_lon, loczoom_lat 592 INTEGER(i_std) :: ndims, nb_atts, id_unlim, lll, iim_full, jjm_full, iim_tmp, jjm_tmp 593 INTEGER(i_std) :: iindex_init, jindex_init, iindex_end, jindex_end, i_orig, j_orig 506 594 CHARACTER(LEN=20) :: varname 507 REAL(r_std),DIMENSION(iim,jjm) :: mapfac_x, mapfac_y595 508 596 INTEGER(i_std), DIMENSION(4) :: vardims 509 597 INTEGER(i_std), DIMENSION(8) :: rose 510 598 ! 599 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_full, lon_full, lat_full, area_full 600 REAL(r_std),ALLOCATABLE, DIMENSION(:,:) :: mapfac_x_full, mapfac_y_full 601 CHARACTER(LEN=20) :: axname 602 REAL(r_std) :: dx, dy, coslat 603 ! 604 REAL(r_std), PARAMETER :: mincos = 0.0001 605 REAL(r_std), PARAMETER :: pi = 3.141592653589793238 606 REAL(r_std), PARAMETER :: R_Earth = 6378000. 607 ! 608 ! A lot of modifications are added to this subroutine (XW) 511 609 ! 512 610 ! Set some default values agains which we can check … … 523 621 calendar = 'gregorian' 524 622 ! 623 ! get dimension of full-grid data 624 ! 625 iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars, & 626 nAttributes=nb_atts, unlimitedDimId=id_unlim) 627 ! 628 DO iv=1,ndims 629 ! 630 iret = NF90_INQUIRE_DIMENSION (fid, iv, name=axname, len=lll) 631 IF (iret /= NF90_NOERR) THEN 632 CALL ipslerr (3,'globgrd_getdomsz',"Could not get size of dimension :"," "," ") 633 ENDIF 634 ! 635 ! This can be refined by testing the actual grid found in the file. 636 ! 637 SELECT CASE(axname) 638 ! 639 !! Coordinate variables used by WRF. 640 CASE("west_east") 641 iim_full = lll 642 643 CASE("south_north") 644 jjm_full = lll 645 646 END SELECT 647 ENDDO 648 525 649 ! 526 650 ! Init projection in grid.f90 so that it can be used later for projections. 527 651 ! 528 CALL grid_initproj(fid, iim , jjm)652 CALL grid_initproj(fid, iim_full, jjm_full) 529 653 ! 530 654 iret = NF90_INQUIRE (fid, nVariables=nvars) … … 533 657 CALL ipslerr (3,'globgrd_getwrf',"Error inquiering variables from WRF grid file."," ", " ") 534 658 ENDIF 659 IF ( .NOT. ALLOCATED(lon_full) ) ALLOCATE(lon_full(iim_full,jjm_full)) 660 IF ( .NOT. ALLOCATED(lat_full) ) ALLOCATE(lat_full(iim_full,jjm_full)) 661 IF ( .NOT. ALLOCATED(mask_full) ) ALLOCATE(mask_full(iim_full,jjm_full)) 662 IF ( .NOT. ALLOCATED(mapfac_x_full) ) ALLOCATE(mapfac_x_full(iim_full,jjm_full)) 663 IF ( .NOT. ALLOCATED(mapfac_y_full) ) ALLOCATE(mapfac_y_full(iim_full,jjm_full)) 535 664 ! 536 665 DO iv = 1,nvars … … 541 670 ! 542 671 CASE("XLONG_M") 543 iret = NF90_GET_VAR(fid, iv, lon )672 iret = NF90_GET_VAR(fid, iv, lon_full) 544 673 IF (iret /= NF90_NOERR) THEN 545 674 CALL ipslerr (3,'globgrd_getwrf',"Could not read the longitude from the WRF grid file.", " ", " ") 546 675 ENDIF 547 676 CASE("XLAT_M") 548 iret = NF90_GET_VAR(fid, iv, lat )677 iret = NF90_GET_VAR(fid, iv, lat_full) 549 678 IF (iret /= NF90_NOERR) THEN 550 679 CALL ipslerr (3,'globgrd_getwrf',"Could not read the latitude from the WRF grid file.", " ", " ") 551 680 ENDIF 552 681 CASE("LANDMASK") 553 iret = NF90_GET_VAR(fid, iv, mask )682 iret = NF90_GET_VAR(fid, iv, mask_full) 554 683 IF (iret /= NF90_NOERR) THEN 555 684 CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ") 556 685 ENDIF 557 686 CASE("MAPFAC_MX") 558 iret = NF90_GET_VAR(fid, iv, mapfac_x )687 iret = NF90_GET_VAR(fid, iv, mapfac_x_full) 559 688 IF (iret /= NF90_NOERR) THEN 560 689 CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ") 561 690 ENDIF 562 691 CASE("MAPFAC_MY") 563 iret = NF90_GET_VAR(fid, iv, mapfac_y )692 iret = NF90_GET_VAR(fid, iv, mapfac_y_full) 564 693 IF (iret /= NF90_NOERR) THEN 565 694 CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ") … … 569 698 ENDDO 570 699 ! 571 ! Compute corners and area on the iimxjjm grid 572 ! 700 IF ( PRESENT(zoom_lon) .AND. PRESENT(zoom_lat) ) THEN 701 ! 702 ! if zoomed region 703 ! 704 loczoom_lon = zoom_lon 705 loczoom_lat = zoom_lat 706 707 708 ! to get new iim, jjm and nbland for the zoomed region 709 ! tmp_lon_ind1, tmp_lat_ind1 etc: dimension(2), 710 ! with the first one for longitude, second one for latitude 711 tmp_lon_ind1 = MINLOC(ABS(lon_full(:,:)-loczoom_lon(1))) 712 tmp_lat_ind1 = MINLOC(ABS(lat_full(:,:)-loczoom_lat(1))) 713 714 tmp_lon_ind2 = MINLOC(ABS(lon_full(:,:)-loczoom_lon(2))) 715 tmp_lat_ind2 = MINLOC(ABS(lat_full(:,:)-loczoom_lat(2))) 716 ! 717 ! get indices for the zoomed region 718 ! lon1, lat2------- lon2,lat2 719 ! | | 720 ! lon1, lat1------- lon2,lat1 721 ! 722 iindex_init = tmp_lon_ind1(1) 723 jindex_init = tmp_lat_ind1(2) 724 725 iindex_end= tmp_lon_ind2(1) 726 jindex_end= tmp_lat_ind2(2) 727 728 ! allocate variable values for the zoomed region 729 mask(:,:) = mask_full(iindex_init:iindex_end, jindex_init:jindex_end) 730 lat(:,:) = lat_full(iindex_init:iindex_end, jindex_init:jindex_end) 731 lon(:,:) = lon_full(iindex_init:iindex_end, jindex_init:jindex_end) 732 733 iim_tmp = iindex_end - iindex_init +1 734 jjm_tmp = jindex_end - jindex_init +1 735 736 737 ELSE 738 ! 739 ! if full grids 740 ! 741 iindex_init = 1 742 jindex_init = 1 743 744 ! define variable values for full grids 745 lat(:,:) = lat_full(:,:) 746 lon(:,:) = lon_full(:,:) 747 mask(:,:) = mask_full(:,:) 748 749 ENDIF 750 ! 751 ! Compute corners on the iimxjjm full or zoomed grid 573 752 DO ip=1,iim 574 753 DO jp=1,jjm 754 i_orig = iindex_init + ip - 1 755 j_orig = jindex_init + jp - 1 575 756 ! Corners 576 CALL grid_tolola(i p+0.5, jp+0.5, corners(ip,jp,1,1), corners(ip,jp,1,2))577 CALL grid_tolola(i p+0.5, jp-0.5, corners(ip,jp,2,1), corners(ip,jp,2,2))578 CALL grid_tolola(i p-0.5, jp-0.5, corners(ip,jp,3,1), corners(ip,jp,3,2))579 CALL grid_tolola(i p-0.5, jp-0.5, corners(ip,jp,4,1), corners(ip,jp,4,2))757 CALL grid_tolola(i_orig+0.5, j_orig+0.5, corners(ip,jp,1,1), corners(ip,jp,1,2)) 758 CALL grid_tolola(i_orig+0.5, j_orig-0.5, corners(ip,jp,2,1), corners(ip,jp,2,2)) 759 CALL grid_tolola(i_orig-0.5, j_orig-0.5, corners(ip,jp,3,1), corners(ip,jp,3,2)) 760 CALL grid_tolola(i_orig-0.5, j_orig-0.5, corners(ip,jp,4,1), corners(ip,jp,4,2)) 580 761 ! 581 762 ENDDO 582 763 ENDDO 583 764 ! 584 ! Compute resolution and area on the gathered grid765 ! Compute resolution and area on the gathered, full or zoomed grid 585 766 ! 586 767 k=0 … … 588 769 DO jp=1,jjm 589 770 DO ip=1,iim 771 ! 772 !get the right index of zoomed/full region in the original grids 773 ! 774 i_orig = iindex_init + ip - 1 775 j_orig = jindex_init + jp - 1 776 ! 777 ! 778 ! compute area 779 coslat = MAX( COS(lat_full(i_orig,j_orig) * pi/180. ), mincos ) 780 dx = ABS(corners(ip,jp,2,1) - corners(ip,jp,1,1)) * pi/180. * R_Earth * coslat 781 dy = ABS(corners(ip,jp,1,2) - corners(ip,jp,3,2)) * pi/180. * R_Earth 782 area(ip,jp) = dx*dy 783 ! 784 ! compute index and contfrac 590 785 IF ( mask(ip,jp) > 0.5 ) THEN 591 786 ! 787 ! index of the points in the local zoomed grid 592 788 k=k+1 593 789 lindex(k) = (jp-1)*iim+ip … … 602 798 CALL ipslerr (3,'globgrd_getwrf',"Error closing the WRF grid file :", " ", " ") 603 799 ENDIF 800 IF ( ALLOCATED(lon_full) ) DEALLOCATE(lon_full) 801 IF ( ALLOCATED(lat_full) ) DEALLOCATE(lat_full) 802 IF ( ALLOCATED(mask_full) ) DEALLOCATE(mask_full) 803 IF ( ALLOCATED(mapfac_x_full) ) DEALLOCATE(mapfac_x_full) 804 IF ( ALLOCATED(mapfac_y_full) ) DEALLOCATE(mapfac_y_full) 805 604 806 ! 605 807 END SUBROUTINE globgrd_getwrf
Note: See TracChangeset
for help on using the changeset viewer.
