Changeset 1275 for trunk/NEMO/OPA_SRC
- Timestamp:
- 2009-01-19T20:22:55+01:00 (15 years ago)
- Location:
- trunk/NEMO/OPA_SRC/SBC
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/SBC/fldread.F90
r1219 r1275 5 5 !!===================================================================== 6 6 !! History : 9.0 ! 06-06 (G. Madec) Original code 7 !! ! 05-08 (S. Alderson) Modified for Interpolation in memory 8 !! ! from input grid to model grid 7 9 !!---------------------------------------------------------------------- 8 10 … … 17 19 USE in_out_manager ! I/O manager 18 20 USE iom ! I/O manager library 21 USE geo2ocean ! for vector rotation on to model grid 19 22 20 23 IMPLICIT NONE … … 28 31 LOGICAL :: ln_clim ! climatology or not (T/F) 29 32 CHARACTER(len = 7) :: cltype ! type of data file 'monthly' or yearly' 33 CHARACTER(len = 34) :: wname ! generic name of a NetCDF weights file to be used, blank if not 34 CHARACTER(len = 34) :: vcomp ! symbolic component name if a vector that needs rotation 35 ! a string starting with "U" or "V" for each component 36 ! chars 2 onwards identify which components go together 30 37 END TYPE FLD_N 31 38 … … 44 51 REAL(wp) , ALLOCATABLE, DIMENSION(:,:) :: fnow ! input fields interpolated to now time step 45 52 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:) :: fdta ! 2 consecutive record of input fields 53 CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file acting as a key 54 ! into the WGTLIST structure 55 CHARACTER(len = 34) :: vcomp ! symbolic name for a vector component that needs rotation 56 LOGICAL , DIMENSION(2) :: rotn ! flag to indicate whether field has been rotated 46 57 END TYPE FLD 58 59 !$AGRIF_DO_NOT_TREAT 60 61 !! keep list of all weights variables so they're only read in once 62 !! need to add AGRIF directives not to process this structure 63 !! also need to force wgtname to include AGRIF nest number 64 TYPE :: WGT !: Input weights related variables 65 CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file 66 INTEGER , DIMENSION(2) :: ddims ! shape of input grid 67 INTEGER , DIMENSION(2) :: botleft ! top left corner of box in input grid containing 68 ! current processor grid 69 INTEGER , DIMENSION(2) :: topright ! top right corner of box 70 INTEGER :: jpiwgt ! width of box on input grid 71 INTEGER :: jpjwgt ! height of box on input grid 72 INTEGER :: numwgt ! number of weights (4=bilinear, 16=bicubic) 73 INTEGER :: nestid ! for agrif, keep track of nest we're in 74 INTEGER :: offset ! =0 when cyclic grid has coincident first/last columns, 75 ! =1 when they assumed to be one grid spacing apart 76 ! =-1 otherwise 77 LOGICAL :: cyclic ! east-west cyclic or not 78 INTEGER, DIMENSION(:,:,:), POINTER :: data_jpi ! array of source integers 79 INTEGER, DIMENSION(:,:,:), POINTER :: data_jpj ! array of source integers 80 REAL(wp), DIMENSION(:,:,:), POINTER :: data_wgt ! array of weights on model grid 81 REAL(wp), DIMENSION(:,:), POINTER :: fly_dta ! array of values on input grid 82 REAL(wp), DIMENSION(:,:), POINTER :: col2 ! temporary array for reading in columns 83 END TYPE WGT 84 85 INTEGER, PARAMETER :: tot_wgts = 10 86 TYPE( WGT ), DIMENSION(tot_wgts) :: ref_wgts ! array of wgts 87 INTEGER :: nxt_wgt = 1 ! point to next available space in ref_wgts array 88 89 !$AGRIF_END_DO_NOT_TREAT 47 90 48 91 PUBLIC fld_read, fld_fill ! called by sbc... modules … … 72 115 TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables 73 116 !! 117 CHARACTER (LEN=34) :: acomp ! dummy weight name 118 INTEGER :: kf, nf ! dummy indices 119 INTEGER :: imf ! size of the structure sd 120 REAL(wp), DIMENSION(jpi,jpj) :: utmp, vtmp! temporary arrays for vector rotation 121 74 122 INTEGER :: jf ! dummy indices 123 INTEGER :: kw ! index into wgts array 75 124 REAL(wp) :: zreclast ! last record to be read in the current year file 76 125 REAL(wp) :: zsecend ! number of second since Jan. 1st 00h of nit000 year at nitend … … 81 130 CHARACTER(LEN=1000) :: clfmt ! write format 82 131 !!--------------------------------------------------------------------- 132 ! 133 imf = SIZE( sd ) 83 134 ! ! ===================== ! 84 DO jf = 1, SIZE( sd )! LOOP OVER FIELD !135 DO jf = 1, imf ! LOOP OVER FIELD ! 85 136 ! ! ===================== ! 86 137 ! … … 93 144 !CDIR COLLAPSE 94 145 sd(jf)%fdta(:,:,1) = sd(jf)%fdta(:,:,2) 146 sd(jf)%rotn(1) = sd(jf)%rotn(2) 95 147 ENDIF 96 148 … … 144 196 145 197 ! read after data 146 CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fdta(:,:,2), NINT( sd(jf)%rec_a(1) ) ) 147 148 ENDIF 149 198 IF( LEN(TRIM(sd(jf)%wgtname)) > 0 ) THEN 199 CALL wgt_list( sd(jf), kw ) 200 CALL fld_interp( sd(jf)%num, sd(jf)%clvar, kw, sd(jf)%fdta(:,:,2), NINT( sd(jf)%rec_a(1) ) ) 201 ELSE 202 CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fdta(:,:,2), NINT( sd(jf)%rec_a(1) ) ) 203 ENDIF 204 sd(jf)%rotn(2) = .FALSE. 205 206 ENDIF 207 ! ! ===================== ! 208 END DO ! END LOOP OVER FIELD ! 209 ! ! ===================== ! 210 211 IF( kt == nit000 .AND. lwp ) CALL wgt_print() 212 213 !! Vector fields may need to be rotated onto the local grid direction 214 !! This has to happen before the time interpolations 215 !! (sga: following code should be modified so that pairs arent searched for each time 216 217 DO jf = 1, imf 218 !! find vector rotations required 219 IF( LEN(TRIM(sd(jf)%vcomp)) > 0 ) THEN 220 !! east-west component has symbolic name starting with 'U' 221 IF( sd(jf)%vcomp(1:1) == 'U' ) THEN 222 !! found an east-west component, look for the north-south component 223 !! which has same symbolic name but with 'U' replaced with 'V' 224 nf = LEN_TRIM( sd(jf)%vcomp ) 225 IF( nf == 1) THEN 226 acomp = 'V' 227 ELSE 228 acomp = 'V' // sd(jf)%vcomp(2:nf) 229 ENDIF 230 kf = -1 231 DO nf = 1, imf 232 IF( TRIM(sd(nf)%vcomp) == TRIM(acomp) ) kf = nf 233 END DO 234 IF( kf > 0 ) THEN 235 !! fields jf,kf are two components which need to be rotated together 236 DO nf = 1,2 237 !! check each time level of this pair 238 IF( .NOT. sd(jf)%rotn(nf) .AND. .NOT. sd(kf)%rotn(nf) ) THEN 239 utmp(:,:) = 0.0 240 vtmp(:,:) = 0.0 241 CALL rot_rep( sd(jf)%fdta(:,:,nf), sd(kf)%fdta(:,:,nf), 'T', 'en->i', utmp(:,:) ) 242 CALL rot_rep( sd(jf)%fdta(:,:,nf), sd(kf)%fdta(:,:,nf), 'T', 'en->j', vtmp(:,:) ) 243 sd(jf)%fdta(:,:,nf) = utmp(:,:) 244 sd(kf)%fdta(:,:,nf) = vtmp(:,:) 245 sd(jf)%rotn(nf) = .TRUE. 246 sd(kf)%rotn(nf) = .TRUE. 247 IF( lwp .AND. kt == nit000 ) & 248 WRITE(numout,*) 'fld_read: vector pair (', & 249 TRIM(sd(jf)%clvar),',',TRIM(sd(kf)%clvar), & 250 ') rotated on to model grid' 251 ENDIF 252 END DO 253 ENDIF 254 ENDIF 255 ENDIF 256 END DO 257 258 ! ! ===================== ! 259 DO jf = 1, imf ! LOOP OVER FIELD ! 260 ! ! ===================== ! 261 ! 150 262 ! update field at each kn_fsbc time-step 151 263 IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN … … 200 312 INTEGER :: idvar ! variable id 201 313 INTEGER :: inrec ! number of record existing for this variable 314 INTEGER :: kwgt 202 315 CHARACTER(LEN=1000) :: clfmt ! write format 203 316 !!--------------------------------------------------------------------- … … 249 362 250 363 ! read before data into sdjf%fdta(:,:,2) because we will swap data in the following part of fld_read 251 CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,2), NINT( sdjf%rec_b(1) ) ) 364 IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 365 CALL wgt_list( sdjf, kwgt ) 366 CALL fld_interp( sdjf%num, sdjf%clvar, kwgt, sdjf%fdta(:,:,2), NINT( sdjf%rec_b(1) ) ) 367 ELSE 368 CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,2), NINT( sdjf%rec_b(1) ) ) 369 ENDIF 370 sdjf%rotn(2) = .FALSE. 252 371 253 372 clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i4, ' at time = ', f7.2, ' days')" … … 418 537 ELSE ; sdf(jf)%cltype = sdf_n(jf)%cltype 419 538 ENDIF 539 sdf(jf)%wgtname = " " 540 IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) & 541 sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname ) 542 sdf(jf)%vcomp = sdf_n(jf)%vcomp 420 543 END DO 421 544 … … 432 555 & ' time interp: ' , sdf(jf)%ln_tint , & 433 556 & ' climatology: ' , sdf(jf)%ln_clim , & 557 & ' weights : ' , TRIM( sdf(jf)%wgtname ), & 558 & ' pairing : ' , TRIM( sdf(jf)%vcomp ), & 434 559 & ' data type: ' , sdf(jf)%cltype 435 560 END DO … … 439 564 440 565 566 SUBROUTINE wgt_list( sd, kwgt ) 567 !!--------------------------------------------------------------------- 568 !! *** ROUTINE wgt_list *** 569 !! 570 !! ** Purpose : search array of WGTs and find a weights file 571 !! entry, or return a new one adding it to the end 572 !! if it is a new entry, the weights data is read in and 573 !! restructured (fld_weight) 574 !! 575 !! ** Method : 576 !!---------------------------------------------------------------------- 577 TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file 578 INTEGER, INTENT(inout) :: kwgt ! index of weights 579 !! 580 INTEGER :: kw 581 INTEGER :: nestid 582 LOGICAL :: found 583 !!---------------------------------------------------------------------- 584 ! 585 !! search down linked list 586 !! weights filename is either present or we hit the end of the list 587 found = .FALSE. 588 589 !! because agrif nest part of filenames are now added in iom_open 590 !! to distinguish between weights files on the different grids, need to track 591 !! nest number explicitly 592 nestid = 0 593 #if defined key_agrif 594 nestid = Agrif_Fixed() 595 #endif 596 DO kw = 1, nxt_wgt-1 597 IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. & 598 ref_wgts(kw)%nestid == nestid) THEN 599 kwgt = kw 600 found = .TRUE. 601 EXIT 602 ENDIF 603 END DO 604 IF( .NOT.found ) THEN 605 kwgt = nxt_wgt 606 CALL fld_weight( sd ) 607 ENDIF 608 609 END SUBROUTINE wgt_list 610 611 SUBROUTINE wgt_print( ) 612 !!--------------------------------------------------------------------- 613 !! *** ROUTINE wgt_print *** 614 !! 615 !! ** Purpose : print the list of known weights 616 !! 617 !! ** Method : 618 !!---------------------------------------------------------------------- 619 !! 620 INTEGER :: kw 621 !!---------------------------------------------------------------------- 622 ! 623 624 DO kw = 1, nxt_wgt-1 625 WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname) 626 WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2) 627 WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt 628 WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt 629 WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt 630 WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft 631 WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright 632 IF( ref_wgts(kw)%cyclic ) THEN 633 WRITE(numout,*) ' cyclical' 634 IF( ref_wgts(kw)%offset > 0 ) WRITE(numout,*) ' with offset' 635 ELSE 636 WRITE(numout,*) ' not cyclical' 637 ENDIF 638 IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated' 639 END DO 640 641 END SUBROUTINE wgt_print 642 643 SUBROUTINE fld_weight( sd ) 644 !!--------------------------------------------------------------------- 645 !! *** ROUTINE fld_weight *** 646 !! 647 !! ** Purpose : create a new WGT structure and fill in data from 648 !! file, restructuring as required 649 !! 650 !! ** Method : 651 !!---------------------------------------------------------------------- 652 TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file 653 !! 654 INTEGER :: jn ! dummy loop indices 655 INTEGER :: inum ! temporary logical unit 656 INTEGER :: id ! temporary variable id 657 CHARACTER (len=5) :: aname 658 INTEGER , DIMENSION(3) :: ddims 659 INTEGER , DIMENSION(jpi, jpj) :: data_src 660 REAL(wp), DIMENSION(jpi, jpj) :: data_tmp 661 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: line2, lines ! temporary array to read 2 lineumns 662 CHARACTER (len=34) :: lonvar 663 LOGICAL :: cyclical 664 REAL(wp) :: resid, dlon ! temporary array to read 2 lineumns 665 INTEGER :: offset ! temporary integer 666 !!---------------------------------------------------------------------- 667 ! 668 IF( nxt_wgt > tot_wgts ) THEN 669 CALL ctl_stop("fld_weights: weights array size exceeded, increase tot_wgts") 670 ENDIF 671 ! 672 !! new weights file entry, add in extra information 673 !! a weights file represents a 2D grid of a certain shape, so we assume that the current 674 !! input data file is representative of all other files to be opened and processed with the 675 !! current weights file 676 677 !! open input data file (non-model grid) 678 CALL iom_open( sd%clname, inum ) 679 680 !! get dimensions 681 id = iom_varid( inum, sd%clvar, ddims ) 682 683 !! check for an east-west cyclic grid 684 !! try to guess name of longitude variable 685 686 lonvar = 'nav_lon' 687 id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.) 688 IF( id <= 0 ) THEN 689 lonvar = 'lon' 690 id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.) 691 ENDIF 692 693 offset = -1 694 cyclical = .FALSE. 695 IF( id > 0 ) THEN 696 !! found a longitude variable 697 !! now going to assume that grid is regular so we can read a single row 698 699 !! because input array is 2d, have to present iom with 2d array even though we only need 1d slice 700 !! worse, we cant pass line2(:,1) to iom_get since this is treated as a 1d array which doesnt match input file 701 ALLOCATE( lines(ddims(1),2) ) 702 CALL iom_get(inum, jpdom_unknown, lonvar, lines(:,:), 1, kstart=(/1,1/), kcount=(/ddims(1),2/) ) 703 704 !! find largest grid spacing 705 lines(1:ddims(1)-1,2) = lines(2:ddims(1),1) - lines(1:ddims(1)-1,1) 706 dlon = MAXVAL( lines(1:ddims(1)-1,2) ) 707 708 resid = ABS(ABS(lines(ddims(1),1)-lines(1,1))-360.0) 709 IF( resid < rsmall ) THEN 710 !! end rows overlap in longitude 711 offset = 0 712 cyclical = .TRUE. 713 ELSEIF( resid < 2.0*dlon ) THEN 714 !! also call it cyclic if difference between end points is less than twice dlon from 360 715 offset = 1 716 cyclical = .TRUE. 717 ENDIF 718 719 DEALLOCATE( lines ) 720 721 ELSE 722 !! guessing failed 723 !! read in first and last columns of data variable 724 !! since we dont know the name of the longitude variable (or even if there is one) 725 !! we assume that if these two columns are equal, file is cyclic east-west 726 727 !! because input array is 2d, have to present iom with 2d array even though we only need 1d slice 728 !! worse, we cant pass line2(1,:) to iom_get since this is treated as a 1d array which doesnt match input file 729 ALLOCATE( lines(2,ddims(2)), line2(2,ddims(2)) ) 730 CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/1,1/), kcount=(/2,ddims(2)/) ) 731 lines(2,:) = line2(1,:) 732 733 CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/ddims(1)-1,1/), kcount=(/2,ddims(2)/) ) 734 lines(1,:) = line2(2,:) 735 736 resid = SUM( ABS(lines(1,:) - lines(2,:)) ) 737 IF( resid < ddims(2)*rsmall ) THEN 738 offset = 0 739 cyclical = .TRUE. 740 ENDIF 741 742 DEALLOCATE( lines, line2 ) 743 ENDIF 744 745 !! close it 746 CALL iom_close( inum ) 747 748 !! now open the weights file 749 750 CALL iom_open ( sd%wgtname, inum ) ! interpolation weights 751 IF ( inum > 0 ) THEN 752 753 ref_wgts(nxt_wgt)%ddims(1) = ddims(1) 754 ref_wgts(nxt_wgt)%ddims(2) = ddims(2) 755 ref_wgts(nxt_wgt)%wgtname = sd%wgtname 756 ref_wgts(nxt_wgt)%offset = -1 757 ref_wgts(nxt_wgt)%cyclic = .FALSE. 758 IF( cyclical ) THEN 759 ref_wgts(nxt_wgt)%offset = offset 760 ref_wgts(nxt_wgt)%cyclic = .TRUE. 761 ENDIF 762 ref_wgts(nxt_wgt)%nestid = 0 763 #if defined key_agrif 764 ref_wgts(nxt_wgt)%nestid = Agrif_Fixed() 765 #endif 766 !! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16) 767 !! for each weight wgtNN there is an integer array srcNN which gives the point in 768 !! the input data grid which is to be multiplied by the weight 769 !! they are both arrays on the model grid so the result of the multiplication is 770 !! added into an output array on the model grid as a running sum 771 772 !! two possible cases: bilinear (4 weights) or bicubic (16 weights) 773 id = iom_varid(inum, 'src05', ldstop=.FALSE.) 774 IF( id <= 0) THEN 775 ref_wgts(nxt_wgt)%numwgt = 4 776 ELSE 777 ref_wgts(nxt_wgt)%numwgt = 16 778 ENDIF 779 780 ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) ) 781 ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) ) 782 ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) ) 783 784 DO jn = 1,4 785 aname = ' ' 786 WRITE(aname,'(a3,i2.2)') 'src',jn 787 data_tmp(:,:) = 0 788 CALL iom_get ( inum, jpdom_unknown, aname, data_tmp(1:nlci,1:nlcj), & 789 kstart=(/nimpp,njmpp/), kcount=(/nlci,nlcj/) ) 790 data_src(:,:) = INT(data_tmp(:,:)) 791 ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1) 792 ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1) 793 END DO 794 795 DO jn = 1, ref_wgts(nxt_wgt)%numwgt 796 aname = ' ' 797 WRITE(aname,'(a3,i2.2)') 'wgt',jn 798 ref_wgts(nxt_wgt)%data_wgt(1:nlci,1:nlcj,jn) = 0.0 799 CALL iom_get ( inum, jpdom_unknown, aname, ref_wgts(nxt_wgt)%data_wgt(1:nlci,1:nlcj,jn), & 800 kstart=(/nimpp,njmpp/), kcount=(/nlci,nlcj/) ) 801 END DO 802 CALL iom_close (inum) 803 804 ! find min and max indices in grid 805 ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(1:nlci,1:nlcj,:)) 806 ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(1:nlci,1:nlcj,:)) 807 ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(1:nlci,1:nlcj,:)) 808 ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(1:nlci,1:nlcj,:)) 809 810 ! and therefore dimensions of the input box 811 ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1 812 ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1 813 814 ! shift indexing of source grid 815 ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1 816 ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1 817 818 ! create input grid, give it a halo to allow gradient calculations 819 ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+2, ref_wgts(nxt_wgt)%jpjwgt+2) ) 820 IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col2(2,ref_wgts(nxt_wgt)%jpjwgt+2) ) 821 822 nxt_wgt = nxt_wgt + 1 823 824 ELSE 825 CALL ctl_stop( ' fld_weight : unable to read the file ' ) 826 ENDIF 827 828 END SUBROUTINE fld_weight 829 830 SUBROUTINE fld_interp(num, clvar, kw, dta, nrec) 831 !!--------------------------------------------------------------------- 832 !! *** ROUTINE fld_interp *** 833 !! 834 !! ** Purpose : apply weights to input gridded data to create data 835 !! on model grid 836 !! 837 !! ** Method : 838 !!---------------------------------------------------------------------- 839 INTEGER, INTENT(in) :: num ! stream number 840 CHARACTER(LEN=*), INTENT(in) :: clvar ! variable name 841 INTEGER, INTENT(in) :: kw ! weights number 842 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: dta ! output field on model grid 843 INTEGER, INTENT(in) :: nrec ! record number to read (ie time slice) 844 !! 845 INTEGER, DIMENSION(2) :: rec1,recn ! temporary arrays for start and length 846 INTEGER :: jk, jn, jm ! loop counters 847 INTEGER :: ni, nj ! lengths 848 INTEGER :: jpimin,jpiwid ! temporary indices 849 INTEGER :: jpjmin,jpjwid ! temporary indices 850 INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices 851 !!---------------------------------------------------------------------- 852 ! 853 854 !! for weighted interpolation we have weights at four corners of a box surrounding 855 !! a model grid point, each weight is multiplied by a grid value (bilinear case) 856 !! or by a grid value and gradients at the corner point (bicubic case) 857 !! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases 858 859 !! sub grid where we already have weights 860 jpimin = ref_wgts(kw)%botleft(1) 861 jpjmin = ref_wgts(kw)%botleft(2) 862 jpiwid = ref_wgts(kw)%jpiwgt 863 jpjwid = ref_wgts(kw)%jpjwgt 864 865 !! what we need to read into sub grid in order to calculate gradients 866 rec1(1) = MAX( jpimin-1, 1 ) 867 rec1(2) = MAX( jpjmin-1, 1 ) 868 recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 ) 869 recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 ) 870 871 !! where we need to read it to 872 jpi1 = 2 + rec1(1) - jpimin 873 jpj1 = 2 + rec1(2) - jpjmin 874 jpi2 = jpi1 + recn(1) - 1 875 jpj2 = jpj1 + recn(2) - 1 876 877 ref_wgts(kw)%fly_dta(:,:) = 0.0 878 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2), nrec, rec1, recn) 879 880 !! first four weights common to both bilinear and bicubic 881 !! note that we have to offset by 1 into fly_dta array because of halo 882 dta(:,:) = 0.0 883 DO jk = 1,4 884 DO jn = 1, jpj 885 DO jm = 1,jpi 886 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 887 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 888 dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1) 889 END DO 890 END DO 891 END DO 892 893 IF (ref_wgts(kw)%numwgt .EQ. 16) THEN 894 895 !! fix up halo points that we couldnt read from file 896 IF( jpi1 == 2 ) THEN 897 ref_wgts(kw)%fly_dta(jpi1-1,:) = ref_wgts(kw)%fly_dta(jpi1,:) 898 ENDIF 899 IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 900 ref_wgts(kw)%fly_dta(jpi2+1,:) = ref_wgts(kw)%fly_dta(jpi2,:) 901 ENDIF 902 IF( jpj1 == 2 ) THEN 903 ref_wgts(kw)%fly_dta(:,jpj1-1) = ref_wgts(kw)%fly_dta(:,jpj1) 904 ENDIF 905 IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN 906 ref_wgts(kw)%fly_dta(:,jpj2+1) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2) - ref_wgts(kw)%fly_dta(:,jpj2-1) 907 ENDIF 908 909 !! if data grid is cyclic we can do better on east-west edges 910 !! but have to allow for whether first and last columns are coincident 911 IF( ref_wgts(kw)%cyclic ) THEN 912 rec1(2) = MAX( jpjmin-1, 1 ) 913 recn(1) = 2 914 recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 ) 915 jpj1 = 2 + rec1(2) - jpjmin 916 jpj2 = jpj1 + recn(2) - 1 917 IF( jpi1 == 2 ) THEN 918 rec1(1) = ref_wgts(kw)%ddims(1) - 1 919 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn) 920 ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2) = ref_wgts(kw)%col2(ref_wgts(kw)%offset+1,jpj1:jpj2) 921 ENDIF 922 IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 923 rec1(1) = 1 924 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn) 925 ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2) = ref_wgts(kw)%col2(2-ref_wgts(kw)%offset,jpj1:jpj2) 926 ENDIF 927 ENDIF 928 929 ! gradient in the i direction 930 DO jk = 1,4 931 DO jn = 1, jpj 932 DO jm = 1,jpi 933 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 934 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 935 dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * & 936 (ref_wgts(kw)%fly_dta(ni+2,nj+1) - ref_wgts(kw)%fly_dta(ni,nj+1)) 937 END DO 938 END DO 939 END DO 940 941 ! gradient in the j direction 942 DO jk = 1,4 943 DO jn = 1, jpj 944 DO jm = 1,jpi 945 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 946 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 947 dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * & 948 (ref_wgts(kw)%fly_dta(ni+1,nj+2) - ref_wgts(kw)%fly_dta(ni+1,nj)) 949 END DO 950 END DO 951 END DO 952 953 ! gradient in the ij direction 954 DO jk = 1,4 955 DO jn = 1, jpj 956 DO jm = 1,jpi 957 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 958 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 959 dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( & 960 (ref_wgts(kw)%fly_dta(ni+2,nj+2) - ref_wgts(kw)%fly_dta(ni ,nj+2)) - & 961 (ref_wgts(kw)%fly_dta(ni+2,nj ) - ref_wgts(kw)%fly_dta(ni ,nj ))) 962 END DO 963 END DO 964 END DO 965 966 END IF 967 968 END SUBROUTINE fld_interp 969 441 970 END MODULE fldread -
trunk/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
r1270 r1275 138 138 139 139 ! (NB: frequency positive => hours, negative => months) 140 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! 141 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! 142 sn_utau = FLD_N( 'utau' , 24. , 'utau' , .true. , .false. , 'yearly' )143 sn_vtau = FLD_N( 'vtau' , 24. , 'vtau' , .true. , .false. , 'yearly' )144 sn_wndm = FLD_N( 'mwnd10m' , 24. , 'm_10' , .true. , .false. , 'yearly' )145 sn_tair = FLD_N( 'tair10m' , 24. , 't_10' , .false. , .false. , 'yearly' )146 sn_humi = FLD_N( 'humi10m' , 24. , 'q_10' , .false. , .false. , 'yearly' )147 sn_ccov = FLD_N( 'ccover' , -1. , 'cloud' , .true. , .false. , 'yearly' )148 sn_prec = FLD_N( 'precip' , -1. , 'precip' , .true. , .false. , ' yearly')140 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! 141 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! 142 sn_utau = FLD_N( 'utau' , 24. , 'utau' , .true. , .false. , 'yearly' , '' , '' ) 143 sn_vtau = FLD_N( 'vtau' , 24. , 'vtau' , .true. , .false. , 'yearly' , '' , '' ) 144 sn_wndm = FLD_N( 'mwnd10m' , 24. , 'm_10' , .true. , .false. , 'yearly' , '' , '' ) 145 sn_tair = FLD_N( 'tair10m' , 24. , 't_10' , .false. , .false. , 'yearly' , '' , '' ) 146 sn_humi = FLD_N( 'humi10m' , 24. , 'q_10' , .false. , .false. , 'yearly' , '' , '' ) 147 sn_ccov = FLD_N( 'ccover' , -1. , 'cloud' , .true. , .false. , 'yearly' , '' , '' ) 148 sn_prec = FLD_N( 'precip' , -1. , 'precip' , .true. , .false. , 'yearly' , '' , '' ) 149 149 150 150 REWIND( numnam ) ! ... read in namlist namsbc_clio -
trunk/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r1270 r1275 125 125 126 126 ! (NB: frequency positive => hours, negative => months) 127 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! 128 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! 129 sn_wndi = FLD_N( 'uwnd10m' , 24. , 'u_10' , .false. , .false. , 'yearly' )130 sn_wndj = FLD_N( 'vwnd10m' , 24. , 'v_10' , .false. , .false. , 'yearly' )131 sn_qsr = FLD_N( 'qsw' , 24. , 'qsw' , .false. , .false. , 'yearly' )132 sn_qlw = FLD_N( 'qlw' , 24. , 'qlw' , .false. , .false. , 'yearly' )133 sn_tair = FLD_N( 'tair10m' , 24. , 't_10' , .false. , .false. , 'yearly' )134 sn_humi = FLD_N( 'humi10m' , 24. , 'q_10' , .false. , .false. , 'yearly' )135 sn_prec = FLD_N( 'precip' , -1. , 'precip' , .true. , .false. , 'yearly' )136 sn_snow = FLD_N( 'snow' , -1. , 'snow' , .true. , .false. , 'yearly' )127 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! 128 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! 129 sn_wndi = FLD_N( 'uwnd10m' , 24. , 'u_10' , .false. , .false. , 'yearly' , '' , '' ) 130 sn_wndj = FLD_N( 'vwnd10m' , 24. , 'v_10' , .false. , .false. , 'yearly' , '' , '' ) 131 sn_qsr = FLD_N( 'qsw' , 24. , 'qsw' , .false. , .false. , 'yearly' , '' , '' ) 132 sn_qlw = FLD_N( 'qlw' , 24. , 'qlw' , .false. , .false. , 'yearly' , '' , '' ) 133 sn_tair = FLD_N( 'tair10m' , 24. , 't_10' , .false. , .false. , 'yearly' , '' , '' ) 134 sn_humi = FLD_N( 'humi10m' , 24. , 'q_10' , .false. , .false. , 'yearly' , '' , '' ) 135 sn_prec = FLD_N( 'precip' , -1. , 'precip' , .true. , .false. , 'yearly' , '' , '' ) 136 sn_snow = FLD_N( 'snow' , -1. , 'snow' , .true. , .false. , 'yearly' , '' , '' ) 137 137 138 138 REWIND( numnam ) ! ... read in namlist namsbc_core -
trunk/NEMO/OPA_SRC/SBC/sbcflx.F90
r1274 r1275 108 108 cn_dir = './' ! directory in which the model is executed 109 109 ! ... default values (NB: frequency positive => hours, negative => months) 110 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! 111 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! 112 sn_utau = FLD_N( 'utau' , 24. , 'utau' , .false. , .false. , 'yearly' )113 sn_vtau = FLD_N( 'vtau' , 24. , 'vtau' , .false. , .false. , 'yearly' )114 sn_qtot = FLD_N( 'qtot' , 24. , 'qtot' , .false. , .false. , 'yearly' )115 sn_qsr = FLD_N( 'qsr' , 24. , 'qsr' , .false. , .false. , 'yearly' )116 sn_emp = FLD_N( 'emp' , 24. , 'emp' , .false. , .false. , 'yearly' )110 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! 111 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! 112 sn_utau = FLD_N( 'utau' , 24. , 'utau' , .false. , .false. , 'yearly' , '' , '' ) 113 sn_vtau = FLD_N( 'vtau' , 24. , 'vtau' , .false. , .false. , 'yearly' , '' , '' ) 114 sn_qtot = FLD_N( 'qtot' , 24. , 'qtot' , .false. , .false. , 'yearly' , '' , '' ) 115 sn_qsr = FLD_N( 'qsr' , 24. , 'qsr' , .false. , .false. , 'yearly' , '' , '' ) 116 sn_emp = FLD_N( 'emp' , 24. , 'emp' , .false. , .false. , 'yearly' , '' , '' ) 117 117 118 118 REWIND ( numnam ) ! ... read in namlist namflx -
trunk/NEMO/OPA_SRC/SBC/sbcice_if.F90
r1200 r1275 69 69 cn_dir = './' ! directory in which the model is executed 70 70 ! ... default values (NB: frequency positive => hours, negative => months) 71 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! 72 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! 73 sn_ice = FLD_N('ice_cover', -1. , 'ice_cov' , .true. , .true. , 'yearly' )71 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! 72 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! 73 sn_ice = FLD_N('ice_cover', -1. , 'ice_cov' , .true. , .true. , 'yearly' , '' , '' ) 74 74 75 75 REWIND ( numnam ) ! ... read in namlist namiif -
trunk/NEMO/OPA_SRC/SBC/sbcrnf.F90
r1242 r1275 128 128 ! ! ============ 129 129 ! (NB: frequency positive => hours, negative => months) 130 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! 131 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! 132 sn_rnf = FLD_N( 'runoffs', -1. , 'sorunoff' , .TRUE. , .true. , 'yearly' )133 sn_cnf = FLD_N( 'runoffs', 0. , 'sorunoff' , .FALSE. , .true. , 'yearly' )130 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! 131 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! 132 sn_rnf = FLD_N( 'runoffs', -1. , 'sorunoff' , .TRUE. , .true. , 'yearly' , '' , '' ) 133 sn_cnf = FLD_N( 'runoffs', 0. , 'sorunoff' , .FALSE. , .true. , 'yearly' , '' , '' ) 134 134 135 135 ! -
trunk/NEMO/OPA_SRC/SBC/sbcssr.F90
r1200 r1275 86 86 cn_dir = './' ! directory in which the model is executed 87 87 ! ... default values (NB: frequency positive => hours, negative => months) 88 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! 89 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! 90 sn_sst = FLD_N( 'sst' , 24. , 'sst' , .false. , .false. , 'yearly' )91 sn_sss = FLD_N( 'sss' , -1. , 'sss' , .true. , .false. , 'yearly' )88 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! 89 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! 90 sn_sst = FLD_N( 'sst' , 24. , 'sst' , .false. , .false. , 'yearly' , '' , '' ) 91 sn_sss = FLD_N( 'sss' , -1. , 'sss' , .true. , .false. , 'yearly' , '' , '' ) 92 92 93 93 REWIND ( numnam ) ! ... read in namlist namflx
Note: See TracChangeset
for help on using the changeset viewer.