- Timestamp:
- 2019-07-09T16:48:48+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11117_obsasm_bugfixes/src/OCE/OBS/obs_grid.F90
r10068 r11230 11 11 !! obs_rlevel_search : Find density level from observed density 12 12 !!---------------------------------------------------------------------- 13 !! * Modules used 13 !! * Modules used 14 14 USE par_kind, ONLY : & ! Precision variables 15 & wp 15 & wp 16 16 USE par_oce, ONLY : & ! Ocean parameters 17 17 & jpk, & … … 32 32 USE netcdf 33 33 USE obs_const, ONLY : & 34 & obfillflt ! Fillvalue 34 & obfillflt ! Fillvalue 35 35 USE lib_mpp, ONLY : & 36 36 & ctl_warn, ctl_stop … … 66 66 REAL(wp), PRIVATE :: dlon ! Lon spacing 67 67 REAL(wp), PRIVATE :: dlat ! Lat spacing 68 68 69 69 INTEGER, PRIVATE :: maxxdiff, maxydiff ! Max diffs between model points 70 70 INTEGER, PRIVATE :: limxdiff, limydiff 71 71 72 72 ! Data storage 73 73 REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: & … … 77 77 & ixpos, & 78 78 & iypos, & 79 & iprocn 79 & iprocn 80 80 81 81 ! Switches … … 83 83 LOGICAL, PUBLIC :: ln_grid_global ! Use global distribution of observations 84 84 CHARACTER(LEN=44), PUBLIC :: & 85 & cn_gridsearchfile ! file name head for grid search lookup 85 & cn_gridsearchfile ! file name head for grid search lookup 86 86 87 87 !!---------------------------------------------------------------------- … … 102 102 !! obs_grd_bruteforce - the original brute force search 103 103 !! or 104 !! obs_grd_lookup - uses a lookup table to do a fast 104 !! obs_grd_lookup - uses a lookup table to do a fast 105 105 !!search 106 106 !!History : 107 !! ! 2007-12 (D. Lea) 107 !! ! 2007-12 (D. Lea) 108 108 !!------------------------------------------------------------------------ 109 109 … … 112 112 & kobsin ! Size of the observation arrays 113 113 REAL(KIND=wp), DIMENSION(kobsin), INTENT(IN) :: & 114 & plam, & ! Longitude of obsrvations 114 & plam, & ! Longitude of obsrvations 115 115 & pphi ! Latitude of observations 116 116 INTEGER, DIMENSION(kobsin), INTENT(OUT) :: & 117 & kobsi, & ! I-index of observations 118 & kobsj, & ! J-index of observations 117 & kobsi, & ! I-index of observations 118 & kobsj, & ! J-index of observations 119 119 & kproc ! Processor number of observations 120 120 CHARACTER(LEN=1) :: & … … 159 159 ENDIF 160 160 ENDIF 161 161 162 162 ENDIF 163 163 … … 165 165 166 166 #include "obs_grd_bruteforce.h90" 167 167 168 168 SUBROUTINE obs_grd_lookup( kobs, plam, pphi, kobsi, kobsj, kproc ) 169 169 !!---------------------------------------------------------------------- … … 177 177 !! ** Action : Return kproc holding the observation and kiobsi,kobsj 178 178 !! valid on kproc=nproc processor only. 179 !! 179 !! 180 180 !! History : 181 181 !! ! 2007-12 (D. Lea) new routine based on obs_grid_search … … 187 187 INTEGER :: kobs ! Size of the observation arrays 188 188 REAL(KIND=wp), DIMENSION(kobs), INTENT(IN) :: & 189 & plam, & ! Longitude of obsrvations 189 & plam, & ! Longitude of obsrvations 190 190 & pphi ! Latitude of observations 191 191 INTEGER, DIMENSION(kobs), INTENT(OUT) :: & 192 & kobsi, & ! I-index of observations 193 & kobsj, & ! J-index of observations 192 & kobsi, & ! I-index of observations 193 & kobsj, & ! J-index of observations 194 194 & kproc ! Processor number of observations 195 195 196 196 !! * Local declarations 197 197 REAL(KIND=wp), DIMENSION(:), ALLOCATABLE :: & … … 301 301 !----------------------------------------------------------------------- 302 302 ! Copy longitudes 303 !----------------------------------------------------------------------- 303 !----------------------------------------------------------------------- 304 304 ALLOCATE( & 305 305 & zplam(kobs) & … … 341 341 DO ji = 1, jlon-1 342 342 zlammax = MAXVAL( zlamtm(:,ji,jj) ) 343 WHERE (zlammax - zlamtm(:, ji, jj) > 180 ) & 343 WHERE (zlammax - zlamtm(:, ji, jj) > 180 ) & 344 344 & zlamtm(:,ji,jj) = zlamtm(:,ji,jj) + 360._wp 345 345 zphitmax(ji,jj) = MAXVAL(zphitm(:,ji,jj)) … … 369 369 ! - For land points kproc is set to number of the processor + 1000000 370 370 ! and we continue the search. 371 ! - For ocean points kproc is set to the number of the processor 371 ! - For ocean points kproc is set to the number of the processor 372 372 ! and we stop the search. 373 373 !----------------------------------------------------------------------- … … 379 379 ! Master loop for grid search 380 380 !------------------------------------------------------------------------ 381 381 382 382 gpkobs: DO jo = 1+joffset, kobs, jostride 383 383 ! Normal case … … 390 390 391 391 ! bottom corner point 392 ipx1 = INT( ( zplam(jo) - lonmin ) / dlon + 1.0 ) 393 ipy1 = INT( ( pphi (jo) - latmin ) / dlat + 1.0 ) 394 392 ipx1 = INT( ( zplam(jo) - lonmin ) / dlon + 1.0 ) 393 ipy1 = INT( ( pphi (jo) - latmin ) / dlat + 1.0 ) 394 395 395 ipx = ipx1 + 1 396 396 ipy = ipy1 + 1 … … 399 399 ! default to false 400 400 llfourflag = .FALSE. 401 401 402 402 ! check for point fully outside of region 403 403 IF ( (ipx1 > nlons) .OR. (ipy1 > nlats) .OR. & … … 415 415 IF (MAXVAL(ixpos(ipx1:ipx,ipy1:ipy)) == -1) CYCLE! cycle if no lookup points found 416 416 ENDIF 417 417 418 418 jimin = 0 419 419 jimax = 0 420 420 jjmin = 0 421 421 jjmax = 0 422 423 IF (.NOT. llfourflag) THEN 422 423 IF (.NOT. llfourflag) THEN 424 424 425 425 ! calculate points range … … 436 436 jjmin = jojjmin - 1 437 437 jjmax = jojjmax + 1 438 439 IF ( jojimin < 0 .OR. jojjmin < 0) THEN 438 439 IF ( jojimin < 0 .OR. jojjmin < 0) THEN 440 440 llfourflag = .TRUE. 441 441 ifourflagcountr(2) = ifourflagcountr(2) + 1 … … 449 449 ifourflagcountr(4) = ifourflagcountr(4) + 1 450 450 ENDIF 451 451 452 452 ENDIF 453 453 … … 455 455 IF (llfourflag) ipmx = 1 456 456 457 IF (llfourflag) THEN 457 IF (llfourflag) THEN 458 458 ifourflagcountt = ifourflagcountt + 1 459 459 ELSE … … 463 463 gridpointsn : DO ip = 0, ipmx 464 464 DO jp = 0, ipmx 465 465 466 466 IF ( kproc(jo) /= -1 ) EXIT gridpointsn 467 467 468 468 ipx = ipx1 + ip 469 469 ipy = ipy1 + jp 470 470 471 471 IF (llfourflag) THEN 472 472 … … 477 477 IF ( ipy < 1 ) ipy = nlats 478 478 479 ! get i,j 479 ! get i,j 480 480 isx = ixpos(ipx,ipy) 481 481 isy = iypos(ipx,ipy) 482 482 483 483 ! estimate appropriate search region (use max/min values) 484 484 jimin = isx - maxxdiff - 1 … … 493 493 IF ( jjmin < 1 ) jjmin = 1 494 494 IF ( jjmax > jlat-1 ) jjmax = jlat-1 495 495 496 496 !--------------------------------------------------------------- 497 ! Ensure that all observation longtiudes are between 0 and 360 497 ! Ensure that all observation longtiudes are between 0 and 360 498 498 !--------------------------------------------------------------- 499 499 500 500 IF ( zplam(jo) < 0.0_wp ) zplam(jo) = zplam(jo) + 360.0_wp 501 501 IF ( zplam(jo) > 360.0_wp ) zplam(jo) = zplam(jo) - 360.0_wp 502 502 503 503 !--------------------------------------------------------------- 504 504 ! Find observations which are on within 1e-6 of a grid point … … 535 535 536 536 IF ( kproc(jo) == -1 ) THEN 537 538 ! Normal case 537 538 ! Normal case 539 539 gridpoints : DO jj = jjmin, jjmax 540 540 DO ji = jimin, jimax … … 543 543 IF ( ( zplam(jo) > zlamtmax(ji,jj) ) .OR. & 544 544 & ( zplam(jo) < zlamtmin(ji,jj) ) ) CYCLE 545 545 546 546 IF ( ABS( pphi(jo) ) < 85 ) THEN 547 547 IF ( ( pphi(jo) > zphitmax(ji,jj) ) .OR. & 548 548 & ( pphi(jo) < zphitmin(ji,jj) ) ) CYCLE 549 549 ENDIF 550 550 551 551 IF ( linquad( zplam(jo), pphi(jo), & 552 552 & zlamtm(:,ji,jj), zphitm(:,ji,jj) ) ) THEN … … 563 563 ENDIF 564 564 ENDIF 565 565 566 566 END DO 567 567 END DO gridpoints … … 572 572 gridpoints_greenwich : DO jj = jjmin, jjmax 573 573 DO ji = jimin, jimax 574 574 575 575 IF ( ( zplam(jo)+360.0_wp > zlamtmax(ji,jj) ) .OR. & 576 576 & ( zplam(jo)+360.0_wp < zlamtmin(ji,jj) ) ) CYCLE … … 595 595 ENDIF 596 596 ENDIF 597 597 598 598 END DO 599 599 END DO gridpoints_greenwich 600 600 601 601 ENDIF ! kproc 602 602 603 603 END DO 604 604 END DO gridpointsn … … 633 633 & zplam & 634 634 & ) 635 635 636 636 END SUBROUTINE obs_grd_lookup 637 637 … … 643 643 !! ** Purpose : Setup a lookup table to reduce the searching required 644 644 !! for converting lat lons to grid point location 645 !! produces or reads in a preexisting file for use in 645 !! produces or reads in a preexisting file for use in 646 646 !! obs_grid_search_lookup_local 647 647 !! 648 !! ** Method : calls obs_grid_search_bruteforce_local with a array 648 !! ** Method : calls obs_grid_search_bruteforce_local with a array 649 649 !! of lats and lons 650 650 !! … … 652 652 !! ! 2007-12 (D. Lea) new routine 653 653 !!---------------------------------------------------------------------- 654 654 655 655 !! * Local declarations 656 656 CHARACTER(LEN=15), PARAMETER :: & 657 657 & cpname = 'obs_grid_setup' 658 CHARACTER(LEN=40) :: cfname 658 CHARACTER(LEN=40) :: cfname 659 659 INTEGER :: ji 660 660 INTEGER :: jj … … 674 674 & lonsi, & 675 675 & latsi 676 INTEGER, DIMENSION(:,:), ALLOCATABLE :: & 676 INTEGER, DIMENSION(:,:), ALLOCATABLE :: & 677 677 & ixposi, & 678 & iyposi, & 679 & iproci 678 & iyposi, & 679 & iproci 680 680 INTEGER, PARAMETER :: histsize=90 681 681 INTEGER, DIMENSION(histsize) :: & … … 684 684 & fhistx1, fhistx2, fhisty1, fhisty2 685 685 REAL(wp) :: histtol 686 686 687 687 IF (ln_grid_search_lookup) THEN 688 688 689 689 WRITE(numout,*) 'Calling obs_grid_setup' 690 690 691 691 IF(lwp) WRITE(numout,*) 692 692 IF(lwp) WRITE(numout,*)'Grid search resolution : ', rn_gridsearchres 693 694 gsearch_nlons_def = NINT( 360.0_wp / rn_gridsearchres ) 693 694 gsearch_nlons_def = NINT( 360.0_wp / rn_gridsearchres ) 695 695 gsearch_nlats_def = NINT( 180.0_wp / rn_gridsearchres ) 696 696 gsearch_lonmin_def = -180.0_wp + 0.5_wp * rn_gridsearchres … … 698 698 gsearch_dlon_def = rn_gridsearchres 699 699 gsearch_dlat_def = rn_gridsearchres 700 700 701 701 IF (lwp) THEN 702 702 WRITE(numout,*)'Grid search gsearch_nlons_def = ',gsearch_nlons_def … … 718 718 fileexist=nf90_open( TRIM( cfname ), nf90_nowrite, & 719 719 & idfile ) 720 720 721 721 IF ( fileexist == nf90_noerr ) THEN 722 722 723 723 ! read data 724 724 ! initially assume size is as defined (to be fixed) 725 725 726 726 WRITE(numout,*) 'Reading: ',cfname 727 727 728 728 CALL chkerr( nf90_open( TRIM( cfname ), nf90_nowrite, idfile ), & 729 729 & cpname, __LINE__ ) 730 730 CALL chkerr( nf90_get_att( idfile, nf90_global, 'maxxdiff', maxxdiff ), & 731 & cpname, __LINE__ ) 731 & cpname, __LINE__ ) 732 732 CALL chkerr( nf90_get_att( idfile, nf90_global, 'maxydiff', maxydiff ), & 733 & cpname, __LINE__ ) 733 & cpname, __LINE__ ) 734 734 CALL chkerr( nf90_get_att( idfile, nf90_global, 'dlon', dlon ), & 735 & cpname, __LINE__ ) 735 & cpname, __LINE__ ) 736 736 CALL chkerr( nf90_get_att( idfile, nf90_global, 'dlat', dlat ), & 737 & cpname, __LINE__ ) 737 & cpname, __LINE__ ) 738 738 CALL chkerr( nf90_get_att( idfile, nf90_global, 'lonmin', lonmin ), & 739 & cpname, __LINE__ ) 739 & cpname, __LINE__ ) 740 740 CALL chkerr( nf90_get_att( idfile, nf90_global, 'latmin', latmin ), & 741 & cpname, __LINE__ ) 742 741 & cpname, __LINE__ ) 742 743 743 CALL chkerr( nf90_inq_dimid(idfile, 'nx' , idnx), & 744 744 & cpname, __LINE__ ) 745 745 CALL chkerr( nf90_inquire_dimension( idfile, idnx, len = nlons ), & 746 & cpname, __LINE__ ) 746 & cpname, __LINE__ ) 747 747 CALL chkerr( nf90_inq_dimid(idfile, 'ny' , idny), & 748 748 & cpname, __LINE__ ) 749 749 CALL chkerr( nf90_inquire_dimension( idfile, idny, len = nlats ), & 750 & cpname, __LINE__ ) 751 750 & cpname, __LINE__ ) 751 752 752 ALLOCATE( & 753 753 & lons(nlons,nlats), & … … 757 757 & iprocn(nlons,nlats) & 758 758 & ) 759 760 CALL chkerr( nf90_inq_varid( idfile, 'XPOS', idxpos ), & 759 760 CALL chkerr( nf90_inq_varid( idfile, 'XPOS', idxpos ), & 761 761 & cpname, __LINE__ ) 762 762 CALL chkerr( nf90_get_var ( idfile, idxpos, ixpos), & 763 763 & cpname, __LINE__ ) 764 CALL chkerr( nf90_inq_varid( idfile, 'YPOS', idypos ), & 764 CALL chkerr( nf90_inq_varid( idfile, 'YPOS', idypos ), & 765 765 & cpname, __LINE__ ) 766 766 CALL chkerr( nf90_get_var ( idfile, idypos, iypos), & 767 767 & cpname, __LINE__ ) 768 768 769 769 CALL chkerr( nf90_close( idfile ), cpname, __LINE__ ) 770 770 771 771 ! setup arrays 772 772 773 773 DO ji = 1, nlons 774 774 DO jj = 1, nlats … … 777 777 END DO 778 778 END DO 779 779 780 780 ! if we are not reading the file we need to create it 781 781 ! create new obs grid search lookup file 782 783 ELSE 784 782 783 ELSE 784 785 785 ! call obs_grid_search 786 786 787 787 IF (lwp) THEN 788 788 WRITE(numout,*) 'creating: ',cfname … … 797 797 dlon = gsearch_dlon_def 798 798 dlat = gsearch_dlat_def 799 799 800 800 ! setup arrays 801 801 802 802 ALLOCATE( & 803 803 & lonsi(nlons,nlats), & … … 807 807 & iproci(nlons,nlats) & 808 808 & ) 809 809 810 810 DO ji = 1, nlons 811 811 DO jj = 1, nlats … … 814 814 END DO 815 815 END DO 816 816 817 817 CALL obs_grd_bruteforce( jpi, jpj, jpiglo, jpjglo, & 818 818 & 1, nlci, 1, nlcj, & … … 821 821 & nlons*nlats, lonsi, latsi, & 822 822 & ixposi, iyposi, iproci ) 823 823 824 824 ! minimise file size by removing regions with no data from xypos file 825 825 ! should be able to just use xpos (ypos will have the same areas of missing data) 826 826 827 827 jimin=1 828 828 jimax=nlons … … 852 852 853 853 maxlat_xpos: DO jj= nlats, 1, -1 854 IF (COUNT(ixposi(:,jj) >= 0) > 0) THEN 854 IF (COUNT(ixposi(:,jj) >= 0) > 0) THEN 855 855 jjmax=jj 856 856 EXIT maxlat_xpos … … 890 890 tmpx2 = 0 891 891 tmpy1 = 0 892 tmpy2 = 0 892 tmpy2 = 0 893 893 894 894 numx1 = 0 … … 933 933 END DO 934 934 935 IF (lwp) THEN936 WRITE(numout,*) 'histograms'937 WRITE(numout,*) '0 1 2 3 4 5 6 7 8 9 10 ...'938 WRITE(numout,*) 'histx1'939 WRITE(numout,*) histx1940 WRITE(numout,*) 'histx2'941 WRITE(numout,*) histx2942 WRITE(numout,*) 'histy1'943 WRITE(numout,*) histy1944 WRITE(numout,*) 'histy2'945 WRITE(numout,*) histy2946 ENDIF947 948 935 meanxdiff1 = tmpx1 / numx1 949 936 meanydiff1 = tmpy1 / numy1 … … 953 940 meanxdiff = MAXVAL((/ meanxdiff1, meanxdiff2 /)) 954 941 meanydiff = MAXVAL((/ meanydiff1, meanydiff2 /)) 955 956 IF (lwp) THEN957 WRITE(numout,*) tmpx1, tmpx2, tmpy1, tmpy2958 WRITE(numout,*) numx1, numx2, numy1, numy2959 WRITE(numout,*) 'meanxdiff: ',meanxdiff, meanxdiff1, meanxdiff2960 WRITE(numout,*) 'meanydiff: ',meanydiff, meanydiff1, meanydiff2961 ENDIF962 942 963 943 tmpx1 = 0 … … 1035 1015 fhisty2(:) = histy2(:) * 1.0 / numy2 1036 1016 1037 ! output new histograms1038 1039 IF (lwp) THEN1040 WRITE(numout,*) 'cumulative histograms'1041 WRITE(numout,*) '0 1 2 3 4 5 6 7 8 9 10 ...'1042 WRITE(numout,*) 'fhistx1'1043 WRITE(numout,*) fhistx11044 WRITE(numout,*) 'fhistx2'1045 WRITE(numout,*) fhistx21046 WRITE(numout,*) 'fhisty1'1047 WRITE(numout,*) fhisty11048 WRITE(numout,*) 'fhisty2'1049 WRITE(numout,*) fhisty21050 ENDIF1051 1052 1017 ! calculate maxxdiff and maxydiff based on cumulative histograms 1053 1018 ! where > 0.999 of points are 1054 1019 1055 ! maxval just converts 1x1 vector return from maxloc to a scalar 1020 ! maxval just converts 1x1 vector return from maxloc to a scalar 1056 1021 1057 1022 histtol = 0.999 … … 1073 1038 CALL chkerr( nf90_put_att( idfile, nf90_global, 'title', & 1074 1039 & 'Mapping file from lon/lat to model grid point' ),& 1075 & cpname,__LINE__ ) 1040 & cpname,__LINE__ ) 1076 1041 CALL chkerr( nf90_put_att( idfile, nf90_global, 'maxxdiff', & 1077 1042 & maxxdiff ), & 1078 & cpname,__LINE__ ) 1043 & cpname,__LINE__ ) 1079 1044 CALL chkerr( nf90_put_att( idfile, nf90_global, 'maxydiff', & 1080 1045 & maxydiff ), & 1081 & cpname,__LINE__ ) 1046 & cpname,__LINE__ ) 1082 1047 CALL chkerr( nf90_put_att( idfile, nf90_global, 'dlon', dlon ),& 1083 & cpname,__LINE__ ) 1048 & cpname,__LINE__ ) 1084 1049 CALL chkerr( nf90_put_att( idfile, nf90_global, 'dlat', dlat ),& 1085 & cpname,__LINE__ ) 1050 & cpname,__LINE__ ) 1086 1051 CALL chkerr( nf90_put_att( idfile, nf90_global, 'lonmin', & 1087 1052 & lonmin ), & 1088 & cpname,__LINE__ ) 1053 & cpname,__LINE__ ) 1089 1054 CALL chkerr( nf90_put_att( idfile, nf90_global, 'latmin', & 1090 1055 & latmin ), & 1091 & cpname,__LINE__ ) 1056 & cpname,__LINE__ ) 1092 1057 1093 1058 CALL chkerr( nf90_def_dim(idfile, 'nx' , nlons, idnx), & … … 1098 1063 incdim(1) = idnx 1099 1064 incdim(2) = idny 1100 1065 1101 1066 CALL chkerr( nf90_def_var( idfile, 'LON', nf90_float, incdim, & 1102 1067 & idlon ), & … … 1105 1070 & 'longitude' ), & 1106 1071 & cpname, __LINE__ ) 1107 1072 1108 1073 CALL chkerr( nf90_def_var( idfile, 'LAT', nf90_float, incdim, & 1109 1074 & idlat ), & … … 1132 1097 1133 1098 CALL chkerr( nf90_enddef( idfile ), cpname, __LINE__ ) 1134 1099 1135 1100 CALL chkerr( nf90_put_var( idfile, idlon, lons), & 1136 1101 & cpname, __LINE__ ) … … 1141 1106 CALL chkerr( nf90_put_var( idfile, idypos, iypos), & 1142 1107 & cpname, __LINE__ ) 1143 1108 1144 1109 CALL chkerr( nf90_close( idfile ), cpname, __LINE__ ) 1145 1146 ! should also output max i, max j spacing for use in 1110 1111 ! should also output max i, max j spacing for use in 1147 1112 ! obs_grid_search_lookup 1148 1113 1149 1114 ENDIF 1150 1115 … … 1154 1119 1155 1120 END SUBROUTINE obs_grid_setup 1156 1121 1157 1122 SUBROUTINE obs_grid_deallocate( ) 1158 1123 !!---------------------------------------------------------------------- … … 1168 1133 DEALLOCATE( lons, lats, ixpos, iypos, iprocn ) 1169 1134 ENDIF 1170 1135 1171 1136 END SUBROUTINE obs_grid_deallocate 1172 1137 … … 1180 1145 1181 1146 END MODULE obs_grid 1182
Note: See TracChangeset
for help on using the changeset viewer.
