Changeset 7646 for trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r6140 r7646 4 4 !! Ocean forcing: read input field for surface boundary condition 5 5 !!===================================================================== 6 !! History : 2.0 ! 06-2006 (S. Masson, G. Madec) Original code 7 !! ! 05-2008 (S. Alderson) Modified for Interpolation in memory from input grid to model grid 8 !! ! 10-2013 (D. Delrosso, P. Oddo) suppression of land point prior to interpolation 6 !! History : 2.0 ! 2006-06 (S. Masson, G. Madec) Original code 7 !! 3.0 ! 2008-05 (S. Alderson) Modified for Interpolation in memory from input grid to model grid 8 !! 3.4 ! 2013-10 (D. Delrosso, P. Oddo) suppression of land point prior to interpolation 9 !! ! 12-2015 (J. Harle) Adding BDY on-the-fly interpolation 9 10 !!---------------------------------------------------------------------- 10 11 11 12 !!---------------------------------------------------------------------- 12 !! fld_read : read input fields used for the computation of the 13 !! surface boundary condition 13 !! fld_read : read input fields used for the computation of the surface boundary condition 14 !! fld_init : initialization of field read 15 !! fld_rec : determined the record(s) to be read 16 !! fld_get : read the data 17 !! fld_map : read global data from file and map onto local data using a general mapping (use for open boundaries) 18 !! fld_rot : rotate the vector fields onto the local grid direction 19 !! fld_clopn : update the data file name and close/open the files 20 !! fld_fill : fill the data structure with the associated information read in namelist 21 !! wgt_list : manage the weights used for interpolation 22 !! wgt_print : print the list of known weights 23 !! fld_weight : create a WGT structure and fill in data from file, restructuring as required 24 !! apply_seaoverland : fill land with ocean values 25 !! seaoverland : create shifted matrices for seaoverland application 26 !! fld_interp : apply weights to input gridded data to create data on model grid 27 !! ksec_week : function returning the first 3 letters of the first day of the weekly file 14 28 !!---------------------------------------------------------------------- 15 29 USE oce ! ocean dynamics and tracers … … 67 81 INTEGER :: nreclast ! last record to be read in the current file 68 82 CHARACTER(len = 256) :: lsmname ! current name of the NetCDF mask file acting as a key 83 INTEGER :: igrd ! grid type for bdy data 84 INTEGER :: ibdy ! bdy set id number 69 85 END TYPE FLD 70 86 … … 114 130 CONTAINS 115 131 116 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset )132 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl ) 117 133 !!--------------------------------------------------------------------- 118 134 !! *** ROUTINE fld_read *** … … 135 151 ! ! kt_offset = +1 => fields at "after" time level 136 152 ! ! etc. 153 INTEGER , INTENT(in ), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 154 LOGICAL , INTENT(in ), OPTIONAL :: fvl ! number of vertical levels in the BDY data 155 !! 137 156 INTEGER :: itmp ! local variable 138 157 INTEGER :: imf ! size of the structure sd … … 171 190 DO jf = 1, imf 172 191 IF( PRESENT(map) ) imap = map(jf) 173 CALL fld_init( kn_fsbc, sd(jf), imap ) ! read each before field (put them in after as they will be swapped) 192 IF( PRESENT(jpk_bdy) ) THEN 193 CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy, fvl ) ! read each before field (put them in after as they will be swapped) 194 ELSE 195 CALL fld_init( kn_fsbc, sd(jf), imap ) ! read each before field (put them in after as they will be swapped) 196 ENDIF 174 197 END DO 175 198 IF( lwp ) CALL wgt_print() ! control print … … 263 286 264 287 ! read after data 265 CALL fld_get( sd(jf), imap ) 266 288 IF( PRESENT(jpk_bdy) ) THEN 289 CALL fld_get( sd(jf), imap, jpk_bdy, fvl) 290 ELSE 291 CALL fld_get( sd(jf), imap ) 292 ENDIF 267 293 ENDIF ! read new data? 268 294 END DO ! --- end loop over field --- ! … … 274 300 IF( sd(jf)%ln_tint ) THEN ! temporal interpolation 275 301 IF(lwp .AND. kt - nit000 <= 100 ) THEN 276 clfmt = "(' fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &302 clfmt = "(' fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // & 277 303 & "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')" 278 304 WRITE(numout, clfmt) TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, & 279 305 & sd(jf)%nrec_b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday 280 WRITE(numout, *) ' it_offset is : ',it_offset306 WRITE(numout, *) ' it_offset is : ',it_offset 281 307 ENDIF 282 308 ! temporal interpolation weights … … 286 312 ELSE ! nothing to do... 287 313 IF(lwp .AND. kt - nit000 <= 100 ) THEN 288 clfmt = "(' fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &314 clfmt = "(' fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // & 289 315 & "', record: ', i6.4, ' (days ', f9.4, ' <-> ', f9.4, ')')" 290 316 WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, & … … 302 328 303 329 304 SUBROUTINE fld_init( kn_fsbc, sdjf, map )330 SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy, fvl) 305 331 !!--------------------------------------------------------------------- 306 332 !! *** ROUTINE fld_init *** … … 309 335 !! - if time interpolation, read before data 310 336 !!---------------------------------------------------------------------- 311 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 312 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 313 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices 337 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 338 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 339 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices 340 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 341 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the BDY data 314 342 !! 315 343 LOGICAL :: llprevyr ! are we reading previous year file? … … 405 433 ENDIF 406 434 ! 407 CALL fld_get( sdjf, map ) ! read before data in after arrays(as we will swap it later) 408 ! 409 clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" 435 ! read before data in after arrays(as we will swap it later) 436 IF( PRESENT(jpk_bdy) ) THEN 437 CALL fld_get( sdjf, map, jpk_bdy, fvl ) 438 ELSE 439 CALL fld_get( sdjf, map ) 440 ENDIF 441 ! 442 clfmt = "(' fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" 410 443 IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday 411 444 ! … … 581 614 582 615 583 SUBROUTINE fld_get( sdjf, map )616 SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl ) 584 617 !!--------------------------------------------------------------------- 585 618 !! *** ROUTINE fld_get *** … … 589 622 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables 590 623 TYPE(MAP_POINTER), INTENT(in ) :: map ! global-to-local mapping indices 624 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the bdy data 625 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the bdy data 591 626 ! 592 627 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) … … 601 636 ! 602 637 IF( ASSOCIATED(map%ptr) ) THEN 603 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map ) 604 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map ) 605 ENDIF 638 IF( PRESENT(jpk_bdy) ) THEN 639 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), & 640 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 641 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), & 642 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 643 ENDIF 644 ELSE 645 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map ) 646 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map ) 647 ENDIF 648 ENDIF 606 649 ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 607 650 CALL wgt_list( sdjf, iw ) … … 658 701 END SUBROUTINE fld_get 659 702 660 661 SUBROUTINE fld_map( num, clvar, dta, nrec, map ) 703 SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl ) 662 704 !!--------------------------------------------------------------------- 663 705 !! *** ROUTINE fld_map *** … … 666 708 !! using a general mapping (for open boundaries) 667 709 !!---------------------------------------------------------------------- 668 #if defined key_bdy 669 USE bdy_oce, ONLY: dta_global, dta_global2! workspace to read in global data arrays670 #endif 710 711 USE bdy_oce, ONLY: ln_bdy, idx_bdy, dta_global, dta_global_z, dta_global_dz, dta_global2, dta_global2_z, dta_global2_dz ! workspace to read in global data arrays 712 671 713 INTEGER , INTENT(in ) :: num ! stream number 672 714 CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name 673 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional)715 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 674 716 INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice) 675 717 TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices 718 INTEGER , INTENT(in), OPTIONAL :: igrd, ibdy, jpk_bdy ! grid type, set number and number of vertical levels in the bdy data 719 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 720 INTEGER :: jpkm1_bdy! number of vertical levels in the bdy data minus 1 676 721 !! 677 722 INTEGER :: ipi ! length of boundary data on local process … … 682 727 INTEGER :: ib, ik, ji, jj ! loop counters 683 728 INTEGER :: ierr 684 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 729 REAL(wp) :: fv ! fillvalue 730 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 731 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read_z ! work space for global data 732 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read_dz ! work space for global data 685 733 !!--------------------------------------------------------------------- 686 734 ! … … 692 740 ilendta = iom_file(num)%dimsz(1,idvar) 693 741 694 #if defined key_bdy 695 ipj = iom_file(num)%dimsz(2,idvar) 696 IF( map%ll_unstruc) THEN ! unstructured open boundary data file 697 dta_read => dta_global 698 ELSE ! structured open boundary data file 699 dta_read => dta_global2 742 IF ( ln_bdy ) THEN 743 ipj = iom_file(num)%dimsz(2,idvar) 744 IF( map%ll_unstruc) THEN ! unstructured open boundary data file 745 dta_read => dta_global 746 IF( PRESENT(jpk_bdy) ) THEN 747 IF( jpk_bdy>0 ) THEN 748 dta_read_z => dta_global_z 749 dta_read_dz => dta_global_dz 750 jpkm1_bdy = jpk_bdy-1 751 ENDIF 752 ENDIF 753 ELSE ! structured open boundary file 754 dta_read => dta_global2 755 IF( PRESENT(jpk_bdy) ) THEN 756 IF( jpk_bdy>0 ) THEN 757 dta_read_z => dta_global2_z 758 dta_read_dz => dta_global2_dz 759 jpkm1_bdy = jpk_bdy-1 760 ENDIF 761 ENDIF 762 ENDIF 700 763 ENDIF 701 #endif702 764 703 765 IF(lwp) WRITE(numout,*) 'Dim size for ', TRIM(clvar),' is ', ilendta … … 705 767 ! 706 768 SELECT CASE( ipk ) 707 CASE(1) ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) 708 CASE DEFAULT ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 769 CASE(1) ; 770 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) 771 IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 772 DO ib = 1, ipi 773 DO ik = 1, ipk 774 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 775 END DO 776 END DO 777 ELSE ! we assume that this is a structured open boundary file 778 DO ib = 1, ipi 779 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 780 ji=map%ptr(ib)-(jj-1)*ilendta 781 DO ik = 1, ipk 782 dta(ib,1,ik) = dta_read(ji,jj,ik) 783 END DO 784 END DO 785 ENDIF 786 787 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 788 ! Do we include something here to adjust barotropic velocities ! 789 ! in case of a depth difference between bdy files and ! 790 ! bathymetry in the case ln_full_vel = .false. and jpk_bdy>0? ! 791 ! [as the enveloping and parital cells could change H] ! 792 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 793 794 CASE DEFAULT ; 795 796 IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN ! boundary data not on model grid: vertical interpolation 797 CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar ) 798 dta_read(:,:,:) = -ABS(fv) 799 dta_read_z(:,:,:) = 0._wp 800 dta_read_dz(:,:,:) = 0._wp 801 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec ) 802 SELECT CASE( igrd ) 803 CASE(1) 804 CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 805 CALL iom_get ( num, jpdom_unknown, 'e3t', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 806 CASE(2) 807 CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 808 CALL iom_get ( num, jpdom_unknown, 'e3u', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 809 CASE(3) 810 CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 811 CALL iom_get ( num, jpdom_unknown, 'e3v', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 812 END SELECT 813 814 IF ( ln_bdy ) & 815 CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 816 817 ELSE ! boundary data assumed to be on model grid 818 819 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 820 IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 821 DO ib = 1, ipi 822 DO ik = 1, ipk 823 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 824 END DO 825 END DO 826 ELSE ! we assume that this is a structured open boundary file 827 DO ib = 1, ipi 828 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 829 ji=map%ptr(ib)-(jj-1)*ilendta 830 DO ik = 1, ipk 831 dta(ib,1,ik) = dta_read(ji,jj,ik) 832 END DO 833 END DO 834 ENDIF 835 ENDIF ! PRESENT(jpk_bdy) 709 836 END SELECT 710 ! 711 IF( map%ll_unstruc ) THEN ! unstructured open boundary data file 837 838 END SUBROUTINE fld_map 839 840 SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 841 842 !!--------------------------------------------------------------------- 843 !! *** ROUTINE fld_bdy_interp *** 844 !! 845 !! ** Purpose : on the fly vertical interpolation to allow the use of 846 !! boundary data from non-native vertical grid 847 !!---------------------------------------------------------------------- 848 USE bdy_oce, ONLY: idx_bdy ! indexing for map <-> ij transformation 849 850 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read ! work space for global data 851 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_z ! work space for global data 852 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_dz ! work space for global data 853 REAL(wp) , INTENT(in) :: fv ! fillvalue and alternative -ABS(fv) 854 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 855 TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices 856 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 857 INTEGER , INTENT(in) :: igrd, ibdy, jpk_bdy ! number of levels in bdy data 858 INTEGER , INTENT(in) :: ilendta ! length of data in file 859 !! 860 INTEGER :: ipi ! length of boundary data on local process 861 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 862 INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 863 INTEGER :: jpkm1_bdy ! number of levels in bdy data minus 1 864 INTEGER :: ib, ik, ikk ! loop counters 865 INTEGER :: ji, jj, zij, zjj ! temporary indices 866 REAL(wp) :: zl, zi, zh ! tmp variable for current depth and interpolation factor 867 REAL(wp) :: fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv) 868 CHARACTER (LEN=10) :: ibstr 869 !!--------------------------------------------------------------------- 870 871 872 ipi = SIZE( dta, 1 ) 873 ipj = SIZE( dta_read, 2 ) 874 ipk = SIZE( dta, 3 ) 875 jpkm1_bdy = jpk_bdy-1 876 877 fv_alt = -ABS(fv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 878 DO ib = 1, ipi 879 zij = idx_bdy(ibdy)%nbi(ib,igrd) 880 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 881 IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj 882 ENDDO 883 ! 884 IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file 885 712 886 DO ib = 1, ipi 713 DO ik = 1, ipk 714 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 887 DO ik = 1, jpk_bdy 888 IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN 889 dta_read_z(map%ptr(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 890 dta_read_dz(map%ptr(ib),1,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 891 ENDIF 892 ENDDO 893 ENDDO 894 895 DO ib = 1, ipi 896 zij = idx_bdy(ibdy)%nbi(ib,igrd) 897 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 898 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 899 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 900 SELECT CASE( igrd ) 901 CASE(1) 902 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 903 WRITE(ibstr,"(I10.10)") map%ptr(ib) 904 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 905 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 906 ENDIF 907 CASE(2) 908 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 909 WRITE(ibstr,"(I10.10)") map%ptr(ib) 910 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 911 IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1), sum(umask(zij,zjj,:)), & 912 & hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1 , & 913 & dta_read(map%ptr(ib),1,:) 914 ENDIF 915 CASE(3) 916 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 917 WRITE(ibstr,"(I10.10)") map%ptr(ib) 918 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 919 ENDIF 920 END SELECT 921 DO ik = 1, ipk 922 SELECT CASE( igrd ) 923 CASE(1) 924 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 925 CASE(2) 926 IF(ln_sco) THEN 927 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 928 ELSE 929 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 930 ENDIF 931 CASE(3) 932 IF(ln_sco) THEN 933 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 934 ELSE 935 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 936 ENDIF 937 END SELECT 938 IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN ! above the first level of external data 939 dta(ib,1,ik) = dta_read(map%ptr(ib),1,1) 940 ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN ! below the last level of external data 941 dta(ib,1,ik) = dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 942 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 943 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 944 IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 945 & .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN 946 zi = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / & 947 & ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) ) 948 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + & 949 & ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi 950 ENDIF 951 END DO 952 ENDIF 715 953 END DO 716 954 END DO 717 ELSE ! structured open boundary data file 955 956 IF(igrd == 2) THEN ! do we need to adjust the transport term? 957 DO ib = 1, ipi 958 zij = idx_bdy(ibdy)%nbi(ib,igrd) 959 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 960 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 961 ztrans = 0._wp 962 ztrans_new = 0._wp 963 DO ik = 1, jpk_bdy ! calculate transport on input grid 964 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 965 ENDDO 966 DO ik = 1, ipk ! calculate transport on model grid 967 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 968 ENDDO 969 DO ik = 1, ipk ! make transport correction 970 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 971 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 972 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 973 IF( ABS(ztrans * r1_hu_n(zij,zjj)) > 0.01_wp ) & 974 & CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 975 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) * umask(zij,zjj,ik) 976 ENDIF 977 ENDDO 978 ENDDO 979 ENDIF 980 981 IF(igrd == 3) THEN ! do we need to adjust the transport term? 982 DO ib = 1, ipi 983 zij = idx_bdy(ibdy)%nbi(ib,igrd) 984 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 985 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 986 ztrans = 0._wp 987 ztrans_new = 0._wp 988 DO ik = 1, jpk_bdy ! calculate transport on input grid 989 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 990 ENDDO 991 DO ik = 1, ipk ! calculate transport on model grid 992 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 993 ENDDO 994 DO ik = 1, ipk ! make transport correction 995 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 996 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 997 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 998 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) * vmask(zij,zjj,ik) 999 ENDIF 1000 ENDDO 1001 ENDDO 1002 ENDIF 1003 1004 ELSE ! structured open boundary file 1005 718 1006 DO ib = 1, ipi 719 1007 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 720 1008 ji=map%ptr(ib)-(jj-1)*ilendta 721 DO ik = 1, ipk 722 dta(ib,1,ik) = dta_read(ji,jj,ik) 1009 DO ik = 1, jpk_bdy 1010 IF( ( dta_read(ji,jj,ik) == fv ) ) THEN 1011 dta_read_z(ji,jj,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 1012 dta_read_dz(ji,jj,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 1013 ENDIF 1014 ENDDO 1015 ENDDO 1016 1017 1018 DO ib = 1, ipi 1019 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1020 ji=map%ptr(ib)-(jj-1)*ilendta 1021 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1022 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1023 zh = SUM(dta_read_dz(ji,jj,:) ) 1024 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 1025 SELECT CASE( igrd ) 1026 CASE(1) 1027 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 1028 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1029 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1030 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 1031 ENDIF 1032 CASE(2) 1033 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 1034 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1035 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1036 ENDIF 1037 CASE(3) 1038 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 1039 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1040 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1041 ENDIF 1042 END SELECT 1043 DO ik = 1, ipk 1044 SELECT CASE( igrd ) ! coded for sco - need zco and zps option using min 1045 CASE(1) 1046 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 1047 CASE(2) 1048 IF(ln_sco) THEN 1049 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1050 ELSE 1051 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 1052 ENDIF 1053 CASE(3) 1054 IF(ln_sco) THEN 1055 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1056 ELSE 1057 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 1058 ENDIF 1059 END SELECT 1060 IF( zl < dta_read_z(ji,jj,1) ) THEN ! above the first level of external data 1061 dta(ib,1,ik) = dta_read(ji,jj,1) 1062 ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN ! below the last level of external data 1063 dta(ib,1,ik) = dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 1064 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 1065 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 1066 IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 1067 & .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN 1068 zi = ( zl - dta_read_z(ji,jj,ikk) ) / & 1069 & ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) ) 1070 dta(ib,1,ik) = dta_read(ji,jj,ikk) + & 1071 & ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi 1072 ENDIF 1073 END DO 1074 ENDIF 723 1075 END DO 724 1076 END DO 725 ENDIF 726 ! 727 END SUBROUTINE fld_map 1077 1078 IF(igrd == 2) THEN ! do we need to adjust the transport term? 1079 DO ib = 1, ipi 1080 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1081 ji=map%ptr(ib)-(jj-1)*ilendta 1082 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1083 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1084 zh = SUM(dta_read_dz(ji,jj,:) ) 1085 ztrans = 0._wp 1086 ztrans_new = 0._wp 1087 DO ik = 1, jpk_bdy ! calculate transport on input grid 1088 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1089 ENDDO 1090 DO ik = 1, ipk ! calculate transport on model grid 1091 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 1092 ENDDO 1093 DO ik = 1, ipk ! make transport correction 1094 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1095 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1096 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1097 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1098 ENDIF 1099 ENDDO 1100 ENDDO 1101 ENDIF 1102 1103 IF(igrd == 3) THEN ! do we need to adjust the transport term? 1104 DO ib = 1, ipi 1105 jj = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1106 ji = map%ptr(ib)-(jj-1)*ilendta 1107 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1108 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1109 zh = SUM(dta_read_dz(ji,jj,:) ) 1110 ztrans = 0._wp 1111 ztrans_new = 0._wp 1112 DO ik = 1, jpk_bdy ! calculate transport on input grid 1113 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1114 ENDDO 1115 DO ik = 1, ipk ! calculate transport on model grid 1116 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 1117 ENDDO 1118 DO ik = 1, ipk ! make transport correction 1119 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1120 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1121 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1122 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1123 ENDIF 1124 ENDDO 1125 ENDDO 1126 ENDIF 1127 1128 ENDIF ! endif unstructured or structured 1129 1130 END SUBROUTINE fld_bdy_interp 728 1131 729 1132 … … 791 1194 !! *** ROUTINE fld_clopn *** 792 1195 !! 793 !! ** Purpose : update the file name and open the file1196 !! ** Purpose : update the file name and close/open the files 794 1197 !!---------------------------------------------------------------------- 795 1198 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables … … 882 1285 883 1286 884 SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )1287 SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam, knoprint ) 885 1288 !!--------------------------------------------------------------------- 886 1289 !! *** ROUTINE fld_fill *** 887 1290 !! 888 !! ** Purpose : fill sdf with sdf_n and control print 1291 !! ** Purpose : fill the data structure (sdf) with the associated information 1292 !! read in namelist (sdf_n) and control print 889 1293 !!---------------------------------------------------------------------- 890 TYPE(FLD) , DIMENSION(:), INTENT(inout) :: sdf ! structure of input fields (file informations, fields read) 891 TYPE(FLD_N), DIMENSION(:), INTENT(in ) :: sdf_n ! array of namelist information structures 892 CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files 893 CHARACTER(len=*) , INTENT(in ) :: cdcaller ! 894 CHARACTER(len=*) , INTENT(in ) :: cdtitle ! 895 CHARACTER(len=*) , INTENT(in ) :: cdnam ! 896 ! 897 INTEGER :: jf ! dummy indices 1294 TYPE(FLD) , DIMENSION(:) , INTENT(inout) :: sdf ! structure of input fields (file informations, fields read) 1295 TYPE(FLD_N), DIMENSION(:) , INTENT(in ) :: sdf_n ! array of namelist information structures 1296 CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files 1297 CHARACTER(len=*) , INTENT(in ) :: cdcaller ! name of the calling routine 1298 CHARACTER(len=*) , INTENT(in ) :: cdtitle ! description of the calling routine 1299 CHARACTER(len=*) , INTENT(in ) :: cdnam ! name of the namelist from which sdf_n comes 1300 INTEGER , OPTIONAL, INTENT(in ) :: knoprint ! no calling routine information printed 1301 ! 1302 INTEGER :: jf ! dummy indices 898 1303 !!--------------------------------------------------------------------- 899 1304 ! … … 922 1327 IF(lwp) THEN ! control print 923 1328 WRITE(numout,*) 924 WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle ) 925 WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /) 926 WRITE(numout,*) ' '//TRIM( cdnam )//' Namelist' 927 WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)' 1329 IF( .NOT.PRESENT( knoprint) ) THEN 1330 WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle ) 1331 WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /) 1332 ENDIF 1333 WRITE(numout,*) ' fld_fill : fill data structure with information from namelist '//TRIM( cdnam ) 1334 WRITE(numout,*) ' ~~~~~~~~' 1335 WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)' 928 1336 DO jf = 1, SIZE(sdf) 929 WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), & 930 & ' variable name: ' , TRIM( sdf(jf)%clvar ) 931 WRITE(numout,*) ' frequency: ' , sdf(jf)%nfreqh , & 932 & ' time interp: ' , sdf(jf)%ln_tint , & 933 & ' climatology: ' , sdf(jf)%ln_clim , & 934 & ' weights : ' , TRIM( sdf(jf)%wgtname ), & 935 & ' pairing : ' , TRIM( sdf(jf)%vcomp ), & 936 & ' data type: ' , sdf(jf)%cltype , & 937 & ' land/sea mask:' , TRIM( sdf(jf)%lsmname ) 1337 WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), ' variable name: ', TRIM( sdf(jf)%clvar ) 1338 WRITE(numout,*) ' frequency: ' , sdf(jf)%nfreqh , & 1339 & ' time interp: ' , sdf(jf)%ln_tint , & 1340 & ' climatology: ' , sdf(jf)%ln_clim 1341 WRITE(numout,*) ' weights: ' , TRIM( sdf(jf)%wgtname ), & 1342 & ' pairing: ' , TRIM( sdf(jf)%vcomp ), & 1343 & ' data type: ' , sdf(jf)%cltype , & 1344 & ' land/sea mask:' , TRIM( sdf(jf)%lsmname ) 938 1345 call flush(numout) 939 1346 END DO … … 947 1354 !! *** ROUTINE wgt_list *** 948 1355 !! 949 !! ** Purpose : search array of WGTs and find a weights file 950 !! entry, or return a new one adding it to the end 951 !! if it is a new entry, the weights data is read in and 952 !! restructured (fld_weight) 1356 !! ** Purpose : search array of WGTs and find a weights file entry, 1357 !! or return a new one adding it to the end if new entry. 1358 !! the weights data is read in and restructured (fld_weight) 953 1359 !!---------------------------------------------------------------------- 954 1360 TYPE( FLD ), INTENT(in ) :: sd ! field with name of weights file … … 1019 1425 !! *** ROUTINE fld_weight *** 1020 1426 !! 1021 !! ** Purpose : create a new WGT structure and fill in data from 1022 !! file,restructuring as required1427 !! ** Purpose : create a new WGT structure and fill in data from file, 1428 !! restructuring as required 1023 1429 !!---------------------------------------------------------------------- 1024 1430 TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file … … 1163 1569 1164 1570 SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm, & 1165 &jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )1571 & jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm ) 1166 1572 !!--------------------------------------------------------------------- 1167 1573 !! *** ROUTINE apply_seaoverland *** … … 1492 1898 !! *** FUNCTION kshift_week *** 1493 1899 !! 1494 !! ** Purpose : 1495 !!--------------------------------------------------------------------- 1496 CHARACTER(len=*), INTENT(in) :: cdday !3 first letters of the first day of the weekly file 1497 !! 1498 INTEGER :: ksec_week ! output variable 1499 INTEGER :: ijul !temp variable 1500 INTEGER :: ishift !temp variable 1900 !! ** Purpose : return the first 3 letters of the first day of the weekly file 1901 !!--------------------------------------------------------------------- 1902 CHARACTER(len=*), INTENT(in) :: cdday ! first 3 letters of the first day of the weekly file 1903 !! 1904 INTEGER :: ksec_week ! output variable 1905 INTEGER :: ijul, ishift ! local integer 1501 1906 CHARACTER(len=3),DIMENSION(7) :: cl_week 1502 1907 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
