Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
- Timestamp:
- 2015-12-21T12:35:23+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r5836 r6140 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 8 !! ! from input grid to model grid 9 !! ! 10-2013 (D. Delrosso, P. Oddo) implement suppression of 10 !! ! land point prior to interpolation 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 11 9 !!---------------------------------------------------------------------- 12 10 … … 15 13 !! surface boundary condition 16 14 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and tracers 18 USE dom_oce ! ocean space and time domain 19 USE phycst ! ??? 20 USE in_out_manager ! I/O manager 21 USE iom ! I/O manager library 22 USE geo2ocean ! for vector rotation on to model grid 23 USE lib_mpp ! MPP library 24 USE wrk_nemo ! work arrays 25 USE lbclnk ! ocean lateral boundary conditions (C1D case) 26 USE ioipsl, ONLY : ymds2ju, ju2ymds ! for calendar 27 USE sbc_oce 15 USE oce ! ocean dynamics and tracers 16 USE dom_oce ! ocean space and time domain 17 USE phycst ! physical constant 18 USE sbc_oce ! surface boundary conditions : fields 19 USE geo2ocean ! for vector rotation on to model grid 20 ! 21 USE in_out_manager ! I/O manager 22 USE iom ! I/O manager library 23 USE ioipsl , ONLY : ymds2ju, ju2ymds ! for calendar 24 USE lib_mpp ! MPP library 25 USE wrk_nemo ! work arrays 26 USE lbclnk ! ocean lateral boundary conditions (C1D case) 28 27 29 28 IMPLICIT NONE … … 60 59 INTEGER , DIMENSION(2) :: nrec_b ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year) 61 60 INTEGER , DIMENSION(2) :: nrec_a ! after record (1: index, 2: second since Jan. 1st 00h of nit000 year) 62 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,: ) :: fnow 63 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) :: fdta 61 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,: ) :: fnow ! input fields interpolated to now time step 62 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) :: fdta ! 2 consecutive record of input fields 64 63 CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file acting as a key 65 64 ! ! into the WGTLIST structure … … 133 132 INTEGER , INTENT(in ), OPTIONAL :: kit ! subcycle timestep for timesplitting option 134 133 INTEGER , INTENT(in ), OPTIONAL :: kt_offset ! provide fields at time other than "now" 135 ! kt_offset = -1 => fields at "before" time level 136 ! kt_offset = +1 => fields at "after" time level 137 ! etc. 138 !! 139 INTEGER :: itmp ! temporary variable 140 INTEGER :: imf ! size of the structure sd 141 INTEGER :: jf ! dummy indices 142 INTEGER :: isecend ! number of second since Jan. 1st 00h of nit000 year at nitend 143 INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step 144 INTEGER :: it_offset ! local time offset variable 145 LOGICAL :: llnxtyr ! open next year file? 146 LOGICAL :: llnxtmth ! open next month file? 147 LOGICAL :: llstop ! stop is the file does not exist 134 ! ! kt_offset = -1 => fields at "before" time level 135 ! ! kt_offset = +1 => fields at "after" time level 136 ! ! etc. 137 INTEGER :: itmp ! local variable 138 INTEGER :: imf ! size of the structure sd 139 INTEGER :: jf ! dummy indices 140 INTEGER :: isecend ! number of second since Jan. 1st 00h of nit000 year at nitend 141 INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step 142 INTEGER :: it_offset ! local time offset variable 143 LOGICAL :: llnxtyr ! open next year file? 144 LOGICAL :: llnxtmth ! open next month file? 145 LOGICAL :: llstop ! stop is the file does not exist 148 146 LOGICAL :: ll_firstcall ! true if this is the first call to fld_read for this set of fields 149 REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation150 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation151 CHARACTER(LEN=1000) :: clfmt 152 TYPE(MAP_POINTER) :: imap ! global-to-local mapping indices147 REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation 148 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation 149 CHARACTER(LEN=1000) :: clfmt ! write format 150 TYPE(MAP_POINTER) :: imap ! global-to-local mapping indices 153 151 !!--------------------------------------------------------------------- 154 152 ll_firstcall = kt == nit000 … … 166 164 isecsbc = nsec_year + nsec1jan000 + (kit+it_offset)*NINT( rdt/REAL(nn_baro,wp) ) 167 165 ELSE ! middle of sbc time step 168 isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdt tra(1)) + it_offset * NINT(rdttra(1))166 isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdt) + it_offset * NINT(rdt) 169 167 ENDIF 170 168 imf = SIZE( sd ) … … 193 191 CALL fld_rec( kn_fsbc, sd(jf), kt_offset = it_offset, kit = kit ) ! update after record informations 194 192 195 ! if kn_fsbc*rdt trais larger than nfreqh (which is kind of odd),193 ! if kn_fsbc*rdt is larger than nfreqh (which is kind of odd), 196 194 ! it is possible that the before value is no more the good one... we have to re-read it 197 195 ! if before is not the last record of the file currently opened and after is the first record to be read … … 214 212 IF( sd(jf)%ln_tint ) THEN 215 213 216 ! if kn_fsbc*rdt trais larger than nfreqh (which is kind of odd),214 ! if kn_fsbc*rdt is larger than nfreqh (which is kind of odd), 217 215 ! it is possible that the before value is no more the good one... we have to re-read it 218 216 ! if before record is not just just before the after record... … … 245 243 ! year/month/week/day file to be not present. If the run continue further than the current 246 244 ! year/month/week/day, next year/month/week/day file must exist 247 isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdt tra(1)) ! second at the end of the run245 isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdt) ! second at the end of the run 248 246 llstop = isecend > sd(jf)%nrec_a(2) ! read more than 1 record of next year 249 247 ! we suppose that the date of next file is next day (should be ok even for weekly files...) … … 299 297 END DO ! --- end loop over field --- ! 300 298 ! 301 ! ! ====================================== ! 302 ENDIF ! update field at each kn_fsbc time-step ! 303 ! ! ====================================== ! 299 ENDIF 304 300 ! 305 301 END SUBROUTINE fld_read … … 333 329 llprevday = .FALSE. 334 330 isec_week = 0 335 331 ! 336 332 ! define record informations 337 333 CALL fld_rec( kn_fsbc, sdjf, ldbefore = .TRUE. ) ! return before values in sdjf%nrec_a (as we will swap it later) 338 334 ! 339 335 ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar 340 336 ! 341 337 IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure 342 338 ! 343 339 IF( sdjf%nrec_a(1) == 0 ) THEN ! we redefine record sdjf%nrec_a(1) with the last record of previous year file 344 340 IF ( sdjf%nfreqh == -12 ) THEN ! yearly mean … … 391 387 ! 392 388 CALL fld_clopn( sdjf, iyear, imonth, iday, .NOT. llprev ) 393 389 ! 394 390 ! if previous year/month/day file does not exist, we switch to the current year/month/day 395 391 IF( llprev .AND. sdjf%num <= 0 ) THEN … … 401 397 CALL fld_clopn( sdjf ) 402 398 ENDIF 403 399 ! 404 400 IF( llprev ) THEN ! check if the record sdjf%nrec_a(1) exists in the file 405 401 idvar = iom_varid( sdjf%num, sdjf%clvar ) ! id of the variable sdjf%clvar … … 408 404 sdjf%nrec_a(1) = MIN( sdjf%nrec_a(1), inrec ) ! make sure we select an existing record 409 405 ENDIF 410 411 ! read before data in after arrays(as we will swap it later) 412 CALL fld_get( sdjf, map ) 413 406 ! 407 CALL fld_get( sdjf, map ) ! read before data in after arrays(as we will swap it later) 408 ! 414 409 clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" 415 410 IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday 416 411 ! 417 412 ENDIF 418 413 ! … … 435 430 LOGICAL , INTENT(in ), OPTIONAL :: ldbefore ! sent back before record values (default = .FALSE.) 436 431 INTEGER , INTENT(in ), OPTIONAL :: kit ! index of barotropic subcycle 437 432 ! ! used only if sdjf%ln_tint = .TRUE. 438 433 INTEGER , INTENT(in ), OPTIONAL :: kt_offset ! Offset of required time level compared to "now" 439 440 ! !434 ! ! time level in units of time steps. 435 ! 441 436 LOGICAL :: llbefore ! local definition of ldbefore 442 437 INTEGER :: iendrec ! end of this record (in seconds) … … 459 454 IF( PRESENT(kt_offset) ) it_offset = kt_offset 460 455 IF( PRESENT(kit) ) THEN ; it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) ) 461 ELSE ; it_offset = it_offset * NINT( rdt tra(1))456 ELSE ; it_offset = it_offset * NINT( rdt ) 462 457 ENDIF 463 458 ! … … 536 531 ELSE ; ztmp = REAL(nsec_year ,wp) ! since 00h on Jan 1 of the current year 537 532 ENDIF 538 ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdt tra(1) + REAL( it_offset, wp )! centrered in the middle of sbc time step539 ztmp = ztmp + 0.01 * rdt tra(1)! avoid truncation error533 ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdt + REAL( it_offset, wp ) ! centrered in the middle of sbc time step 534 ztmp = ztmp + 0.01 * rdt ! avoid truncation error 540 535 IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record 541 536 ! … … 592 587 !! ** Purpose : read the data 593 588 !!---------------------------------------------------------------------- 594 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables595 TYPE(MAP_POINTER), INTENT(in) :: map! global-to-local mapping indices596 ! !597 INTEGER :: ipk! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )598 INTEGER :: iw! index into wgts array599 INTEGER :: ipdom! index of the domain600 INTEGER :: idvar! variable ID601 INTEGER :: idmspc! number of spatial dimensions602 LOGICAL :: lmoor! C1D case: point data589 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables 590 TYPE(MAP_POINTER), INTENT(in ) :: map ! global-to-local mapping indices 591 ! 592 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 593 INTEGER :: iw ! index into wgts array 594 INTEGER :: ipdom ! index of the domain 595 INTEGER :: idvar ! variable ID 596 INTEGER :: idmspc ! number of spatial dimensions 597 LOGICAL :: lmoor ! C1D case: point data 603 598 !!--------------------------------------------------------------------- 604 599 ! … … 611 606 ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 612 607 CALL wgt_list( sdjf, iw ) 613 IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fdta(:,:,:,2),&614 &sdjf%nrec_a(1), sdjf%lsmname )615 ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fnow(:,:,: ),&616 &sdjf%nrec_a(1), sdjf%lsmname )608 IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fdta(:,:,:,2), & 609 & sdjf%nrec_a(1), sdjf%lsmname ) 610 ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fnow(:,:,: ), & 611 & sdjf%nrec_a(1), sdjf%lsmname ) 617 612 ENDIF 618 613 ELSE 619 IF( SIZE(sdjf%fnow, 1) == jpi ) THEN ;ipdom = jpdom_data620 ELSE ;ipdom = jpdom_unknown614 IF( SIZE(sdjf%fnow, 1) == jpi ) THEN ; ipdom = jpdom_data 615 ELSE ; ipdom = jpdom_unknown 621 616 ENDIF 622 617 ! C1D case: If product of spatial dimensions == ipk, then x,y are of 623 618 ! size 1 (point/mooring data): this must be read onto the central grid point 624 619 idvar = iom_varid( sdjf%num, sdjf%clvar ) 625 idmspc = iom_file ( sdjf%num )%ndims( idvar )620 idmspc = iom_file ( sdjf%num )%ndims( idvar ) 626 621 IF( iom_file( sdjf%num )%luld( idvar ) ) idmspc = idmspc - 1 627 lmoor = ( idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk)622 lmoor = ( idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk ) 628 623 ! 629 624 SELECT CASE( ipk ) … … 660 655 ! 661 656 sdjf%rotn(2) = .false. ! vector not yet rotated 662 657 ! 663 658 END SUBROUTINE fld_get 659 664 660 665 661 SUBROUTINE fld_map( num, clvar, dta, nrec, map ) … … 688 684 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 689 685 !!--------------------------------------------------------------------- 690 686 ! 691 687 ipi = SIZE( dta, 1 ) 692 688 ipj = 1 693 689 ipk = SIZE( dta, 3 ) 694 690 ! 695 691 idvar = iom_varid( num, clvar ) 696 692 ilendta = iom_file(num)%dimsz(1,idvar) … … 698 694 #if defined key_bdy 699 695 ipj = iom_file(num)%dimsz(2,idvar) 700 IF ( map%ll_unstruc) THEN! unstructured open boundary data file696 IF( map%ll_unstruc) THEN ! unstructured open boundary data file 701 697 dta_read => dta_global 702 ELSE ! structured open boundary data file698 ELSE ! structured open boundary data file 703 699 dta_read => dta_global2 704 700 ENDIF 705 701 #endif 706 702 707 IF(lwp) WRITE(numout,*) 'Dim size for ', TRIM(clvar),' is ', ilendta703 IF(lwp) WRITE(numout,*) 'Dim size for ', TRIM(clvar),' is ', ilendta 708 704 IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk 709 705 ! 710 706 SELECT CASE( ipk ) 711 707 CASE(1) ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) … … 713 709 END SELECT 714 710 ! 715 IF ( map%ll_unstruc ) THEN! unstructured open boundary data file711 IF( map%ll_unstruc ) THEN ! unstructured open boundary data file 716 712 DO ib = 1, ipi 717 713 DO ik = 1, ipk … … 728 724 END DO 729 725 ENDIF 730 726 ! 731 727 END SUBROUTINE fld_map 732 728 … … 738 734 !! ** Purpose : Vector fields may need to be rotated onto the local grid direction 739 735 !!---------------------------------------------------------------------- 740 INTEGER , INTENT(in ) :: kt ! ocean time step 741 TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables 742 !! 743 INTEGER :: ju,jv,jk,jn ! loop indices 744 INTEGER :: imf ! size of the structure sd 745 INTEGER :: ill ! character length 746 INTEGER :: iv ! indice of V component 736 INTEGER , INTENT(in ) :: kt ! ocean time step 737 TYPE(FLD), DIMENSION(:), INTENT(inout) :: sd ! input field related variables 738 ! 739 INTEGER :: ju, jv, jk, jn ! loop indices 740 INTEGER :: imf ! size of the structure sd 741 INTEGER :: ill ! character length 742 INTEGER :: iv ! indice of V component 743 CHARACTER (LEN=100) :: clcomp ! dummy weight name 747 744 REAL(wp), POINTER, DIMENSION(:,:) :: utmp, vtmp ! temporary arrays for vector rotation 748 CHARACTER (LEN=100) :: clcomp ! dummy weight name 749 !!--------------------------------------------------------------------- 750 751 CALL wrk_alloc( jpi,jpj, utmp, vtmp ) 752 745 !!--------------------------------------------------------------------- 746 ! 747 CALL wrk_alloc( jpi,jpj, utmp, vtmp ) 748 ! 753 749 !! (sga: following code should be modified so that pairs arent searched for each time 754 750 ! … … 786 782 END DO 787 783 ! 788 CALL wrk_dealloc( jpi,jpj, utmp, vtmp )784 CALL wrk_dealloc( jpi,jpj, utmp, vtmp ) 789 785 ! 790 786 END SUBROUTINE fld_rot … … 802 798 INTEGER, OPTIONAL, INTENT(in ) :: kday ! day value 803 799 LOGICAL, OPTIONAL, INTENT(in ) :: ldstop ! stop if open to read a non-existing file (default = .TRUE.) 804 ! !800 ! 805 801 LOGICAL :: llprevyr ! are we reading previous year file? 806 802 LOGICAL :: llprevmth ! are we reading previous month file? … … 853 849 ! 854 850 IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN ! new file to be open 855 851 ! 856 852 sdjf%clname = TRIM(clname) 857 853 IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open 858 854 CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 ) 859 855 ! 860 856 ! find the last record to be read -> update sdjf%nreclast 861 857 indexyr = iyear - nyear + 1 … … 866 862 CASE ( 2 ) ; imonth_len = 31 ! next year -> imonth = 1 867 863 END SELECT 868 864 ! 869 865 ! last record to be read in the current file 870 866 IF ( sdjf%nfreqh == -12 ) THEN ; sdjf%nreclast = 1 ! yearly mean … … 880 876 ENDIF 881 877 ENDIF 882 878 ! 883 879 ENDIF 884 880 ! … … 901 897 INTEGER :: jf ! dummy indices 902 898 !!--------------------------------------------------------------------- 903 899 ! 904 900 DO jf = 1, SIZE(sdf) 905 901 sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname ) … … 923 919 sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 924 920 END DO 925 921 ! 926 922 IF(lwp) THEN ! control print 927 923 WRITE(numout,*) … … 943 939 END DO 944 940 ENDIF 945 941 ! 946 942 END SUBROUTINE fld_fill 947 943 … … 958 954 TYPE( FLD ), INTENT(in ) :: sd ! field with name of weights file 959 955 INTEGER , INTENT(inout) :: kwgt ! index of weights 960 ! !956 ! 961 957 INTEGER :: kw, nestid ! local integer 962 958 LOGICAL :: found ! local logical … … 966 962 !! weights filename is either present or we hit the end of the list 967 963 found = .FALSE. 968 964 ! 969 965 !! because agrif nest part of filenames are now added in iom_open 970 966 !! to distinguish between weights files on the different grids, need to track … … 1028 1024 TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file 1029 1025 !! 1030 INTEGER :: jn ! dummy loop indices 1031 INTEGER :: inum ! temporary logical unit 1032 INTEGER :: id ! temporary variable id 1033 INTEGER :: ipk ! temporary vertical dimension 1034 CHARACTER (len=5) :: aname 1026 INTEGER :: jn ! dummy loop indices 1027 INTEGER :: inum ! local logical unit 1028 INTEGER :: id ! local variable id 1029 INTEGER :: ipk ! local vertical dimension 1030 INTEGER :: zwrap ! local integer 1031 LOGICAL :: cyclical ! 1032 CHARACTER (len=5) :: aname ! 1035 1033 INTEGER , DIMENSION(:), ALLOCATABLE :: ddims 1036 1034 INTEGER , POINTER, DIMENSION(:,:) :: data_src 1037 1035 REAL(wp), POINTER, DIMENSION(:,:) :: data_tmp 1038 LOGICAL :: cyclical 1039 INTEGER :: zwrap ! local integer 1040 !!---------------------------------------------------------------------- 1041 ! 1042 CALL wrk_alloc( jpi,jpj, data_src ) ! integer 1043 CALL wrk_alloc( jpi,jpj, data_tmp ) 1036 !!---------------------------------------------------------------------- 1037 ! 1038 CALL wrk_alloc( jpi,jpj, data_src ) ! integer 1039 CALL wrk_alloc( jpi,jpj, data_tmp ) 1044 1040 ! 1045 1041 IF( nxt_wgt > tot_wgts ) THEN … … 1151 1147 ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) ) 1152 1148 IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) ) 1153 1149 ! 1154 1150 nxt_wgt = nxt_wgt + 1 1155 1151 ! 1156 1152 ELSE 1157 1153 CALL ctl_stop( ' fld_weight : unable to read the file ' ) … … 1166 1162 1167 1163 1168 SUBROUTINE apply_seaoverland( clmaskfile,zfieldo,jpi1_lsm,jpi2_lsm,jpj1_lsm,&1169 & jpj2_lsm, itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)1164 SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm, & 1165 & jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm ) 1170 1166 !!--------------------------------------------------------------------- 1171 1167 !! *** ROUTINE apply_seaoverland *** … … 1176 1172 !! D. Delrosso INGV 1177 1173 !!---------------------------------------------------------------------- 1178 INTEGER :: inum,jni,jnj,jnz,jc ! temporary indices1179 INTEGER, INTENT(in ) :: itmpi,itmpj,itmpz ! lengths1180 INTEGER, INTENT(in) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices1181 INTEGER, DIMENSION(3), INTENT(in) :: rec1_lsm,recn_lsm ! temporary arrays for start and length1182 REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo ! input/output array for seaoverland application1183 REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1 ! temporary array for land point detection1184 REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfieldn ! array of forcing field with undeff for land points1185 REAL(wp),DIMENSION (:,: ), ALLOCATABLE :: zfield ! array of forcing field1186 CHARACTER (len=100), INTENT(in) :: clmaskfile ! land/sea mask file name1187 !!---------------------------------------------------------------------1188 ALLOCATE ( zslmec1(itmpi,itmpj,itmpz) )1189 ALLOCATE ( zfieldn(itmpi,itmpj) )1190 ALLOCATE ( z field(itmpi,itmpj) )1191 1174 INTEGER, INTENT(in ) :: itmpi,itmpj,itmpz ! lengths 1175 INTEGER, INTENT(in ) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices 1176 INTEGER, DIMENSION(3), INTENT(in ) :: rec1_lsm,recn_lsm ! temporary arrays for start and length 1177 REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo ! input/output array for seaoverland application 1178 CHARACTER (len=100), INTENT(in ) :: clmaskfile ! land/sea mask file name 1179 ! 1180 INTEGER :: inum,jni,jnj,jnz,jc ! local indices 1181 REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1 ! local array for land point detection 1182 REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfieldn ! array of forcing field with undeff for land points 1183 REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfield ! array of forcing field 1184 !!--------------------------------------------------------------------- 1185 ! 1186 ALLOCATE ( zslmec1(itmpi,itmpj,itmpz), zfieldn(itmpi,itmpj), zfield(itmpi,itmpj) ) 1187 ! 1192 1188 ! Retrieve the land sea mask data 1193 1189 CALL iom_open( clmaskfile, inum ) 1194 1190 SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) ) 1195 1191 CASE(1) 1196 1192 CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm) 1197 1193 CASE DEFAULT 1198 1194 CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm) 1199 1195 END SELECT 1200 1196 CALL iom_close( inum ) 1201 1202 DO jnz=1,rec1_lsm(3) !! Loop over k dimension 1203 1204 DO jni=1,itmpi !! copy the original field into a tmp array 1205 DO jnj=1,itmpj !! substituting undeff over land points 1206 zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz) 1207 IF ( zslmec1(jni,jnj,jnz) == 1. ) THEN 1208 zfieldn(jni,jnj) = undeff_lsm 1209 ENDIF 1197 ! 1198 DO jnz=1,rec1_lsm(3) !! Loop over k dimension 1199 ! 1200 DO jni = 1, itmpi !! copy the original field into a tmp array 1201 DO jnj = 1, itmpj !! substituting undeff over land points 1202 zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz) 1203 IF( zslmec1(jni,jnj,jnz) == 1. ) zfieldn(jni,jnj) = undeff_lsm 1210 1204 END DO 1211 1205 END DO 1212 1213 CALL seaoverland(zfieldn,itmpi,itmpj,zfield) 1214 DO jc=1,nn_lsm 1215 CALL seaoverland(zfield,itmpi,itmpj,zfield) 1216 END DO 1217 1218 ! Check for Undeff and substitute original values 1219 IF(ANY(zfield==undeff_lsm)) THEN 1220 DO jni=1,itmpi 1221 DO jnj=1,itmpj 1222 IF (zfield(jni,jnj)==undeff_lsm) THEN 1223 zfield(jni,jnj) = zfieldo(jni,jnj,jnz) 1224 ENDIF 1225 ENDDO 1226 ENDDO 1227 ENDIF 1228 1229 zfieldo(:,:,jnz)=zfield(:,:) 1230 1231 END DO !! End Loop over k dimension 1232 1233 DEALLOCATE ( zslmec1 ) 1234 DEALLOCATE ( zfieldn ) 1235 DEALLOCATE ( zfield ) 1236 1206 ! 1207 CALL seaoverland( zfieldn, itmpi, itmpj, zfield ) 1208 DO jc = 1, nn_lsm 1209 CALL seaoverland( zfield, itmpi, itmpj, zfield ) 1210 END DO 1211 ! 1212 ! Check for Undeff and substitute original values 1213 IF( ANY(zfield==undeff_lsm) ) THEN 1214 DO jni = 1, itmpi 1215 DO jnj = 1, itmpj 1216 IF( zfield(jni,jnj)==undeff_lsm ) zfield(jni,jnj) = zfieldo(jni,jnj,jnz) 1217 END DO 1218 END DO 1219 ENDIF 1220 ! 1221 zfieldo(:,:,jnz) = zfield(:,:) 1222 ! 1223 END DO !! End Loop over k dimension 1224 ! 1225 DEALLOCATE ( zslmec1, zfieldn, zfield ) 1226 ! 1237 1227 END SUBROUTINE apply_seaoverland 1238 1228 1239 1229 1240 SUBROUTINE seaoverland( zfieldn,ileni,ilenj,zfield)1230 SUBROUTINE seaoverland( zfieldn, ileni, ilenj, zfield ) 1241 1231 !!--------------------------------------------------------------------- 1242 1232 !! *** ROUTINE seaoverland *** … … 1245 1235 !! D. Delrosso INGV 1246 1236 !!---------------------------------------------------------------------- 1247 INTEGER,INTENT(in) :: ileni,ilenj ! lengths 1248 REAL,DIMENSION (ileni,ilenj),INTENT(in) :: zfieldn ! array of forcing field with undeff for land points 1249 REAL,DIMENSION (ileni,ilenj),INTENT(out) :: zfield ! array of forcing field 1250 REAL,DIMENSION (ileni,ilenj) :: zmat1,zmat2,zmat3,zmat4 ! temporary arrays for seaoverland application 1251 REAL,DIMENSION (ileni,ilenj) :: zmat5,zmat6,zmat7,zmat8 ! temporary arrays for seaoverland application 1252 REAL,DIMENSION (ileni,ilenj) :: zlsm2d ! temporary arrays for seaoverland application 1253 REAL,DIMENSION (ileni,ilenj,8) :: zlsm3d ! temporary arrays for seaoverland application 1254 LOGICAL,DIMENSION (ileni,ilenj,8) :: ll_msknan3d ! logical mask for undeff detection 1255 LOGICAL,DIMENSION (ileni,ilenj) :: ll_msknan2d ! logical mask for undeff detection 1237 INTEGER , INTENT(in ) :: ileni,ilenj ! lengths 1238 REAL, DIMENSION (ileni,ilenj), INTENT(in ) :: zfieldn ! array of forcing field with undeff for land points 1239 REAL, DIMENSION (ileni,ilenj), INTENT( out) :: zfield ! array of forcing field 1240 ! 1241 REAL , DIMENSION (ileni,ilenj) :: zmat1, zmat2, zmat3, zmat4 ! local arrays 1242 REAL , DIMENSION (ileni,ilenj) :: zmat5, zmat6, zmat7, zmat8 ! - - 1243 REAL , DIMENSION (ileni,ilenj) :: zlsm2d ! - - 1244 REAL , DIMENSION (ileni,ilenj,8) :: zlsm3d ! - - 1245 LOGICAL, DIMENSION (ileni,ilenj,8) :: ll_msknan3d ! logical mask for undeff detection 1246 LOGICAL, DIMENSION (ileni,ilenj) :: ll_msknan2d ! logical mask for undeff detection 1256 1247 !!---------------------------------------------------------------------- 1257 zmat8 = eoshift( zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(:,1)/) ,DIM=2)1258 zmat1 = eoshift( zmat8 , SHIFT=-1, BOUNDARY = (/zmat8(1,:)/) ,DIM=1)1259 zmat2 = eoshift( zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(1,:)/) ,DIM=1)1260 zmat4 = eoshift( zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(:,ilenj)/),DIM=2)1261 zmat3 = eoshift( zmat4 , SHIFT=-1, BOUNDARY = (/zmat4(1,:)/) ,DIM=1)1262 zmat5 = eoshift( zmat4 , SHIFT= 1, BOUNDARY = (/zmat4(ileni,:)/) ,DIM=1)1263 zmat6 = eoshift( zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(ileni,:)/),DIM=1)1264 zmat7 = eoshift( zmat8 , SHIFT= 1, BOUNDARY = (/zmat8(ileni,:)/) ,DIM=1)1265 1248 zmat8 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(:,1)/) , DIM=2 ) 1249 zmat1 = eoshift( zmat8 , SHIFT=-1 , BOUNDARY = (/zmat8(1,:)/) , DIM=1 ) 1250 zmat2 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(1,:)/) , DIM=1 ) 1251 zmat4 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(:,ilenj)/) , DIM=2 ) 1252 zmat3 = eoshift( zmat4 , SHIFT=-1 , BOUNDARY = (/zmat4(1,:)/) , DIM=1 ) 1253 zmat5 = eoshift( zmat4 , SHIFT= 1 , BOUNDARY = (/zmat4(ileni,:)/) , DIM=1 ) 1254 zmat6 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(ileni,:)/) , DIM=1 ) 1255 zmat7 = eoshift( zmat8 , SHIFT= 1 , BOUNDARY = (/zmat8(ileni,:)/) , DIM=1 ) 1256 ! 1266 1257 zlsm3d = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /)) 1267 ll_msknan3d = .not.(zlsm3d==undeff_lsm) 1268 ll_msknan2d = .not.(zfieldn==undeff_lsm) ! FALSE where is Undeff (land) 1269 zlsm2d = (SUM ( zlsm3d, 3 , ll_msknan3d ) )/(MAX(1,(COUNT( ll_msknan3d , 3 )) )) 1270 WHERE ((COUNT( ll_msknan3d , 3 )) == 0.0_wp) zlsm2d = undeff_lsm 1271 zfield = MERGE (zfieldn,zlsm2d,ll_msknan2d) 1258 ll_msknan3d = .NOT.( zlsm3d == undeff_lsm ) 1259 ll_msknan2d = .NOT.( zfieldn == undeff_lsm ) ! FALSE where is Undeff (land) 1260 zlsm2d = SUM( zlsm3d, 3 , ll_msknan3d ) / MAX( 1 , COUNT( ll_msknan3d , 3 ) ) 1261 WHERE( COUNT( ll_msknan3d , 3 ) == 0._wp ) zlsm2d = undeff_lsm 1262 zfield = MERGE( zfieldn, zlsm2d, ll_msknan2d ) 1263 ! 1272 1264 END SUBROUTINE seaoverland 1273 1265 … … 1288 1280 INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice) 1289 1281 CHARACTER(LEN=*) , INTENT(in ) :: lsmfile ! land sea mask file name 1290 !! 1291 REAL(wp),DIMENSION(:,:,:),ALLOCATABLE :: ztmp_fly_dta,zfieldo ! temporary array of values on input grid 1292 INTEGER, DIMENSION(3) :: rec1,recn ! temporary arrays for start and length 1293 INTEGER, DIMENSION(3) :: rec1_lsm,recn_lsm ! temporary arrays for start and length in case of seaoverland 1294 INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices 1295 INTEGER :: jk, jn, jm, jir, jjr ! loop counters 1296 INTEGER :: ni, nj ! lengths 1297 INTEGER :: jpimin,jpiwid ! temporary indices 1298 INTEGER :: jpimin_lsm,jpiwid_lsm ! temporary indices 1299 INTEGER :: jpjmin,jpjwid ! temporary indices 1300 INTEGER :: jpjmin_lsm,jpjwid_lsm ! temporary indices 1301 INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices 1302 INTEGER :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices 1303 INTEGER :: itmpi,itmpj,itmpz ! lengths 1304 1282 ! 1283 INTEGER, DIMENSION(3) :: rec1, recn ! temporary arrays for start and length 1284 INTEGER, DIMENSION(3) :: rec1_lsm, recn_lsm ! temporary arrays for start and length in case of seaoverland 1285 INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices 1286 INTEGER :: jk, jn, jm, jir, jjr ! loop counters 1287 INTEGER :: ni, nj ! lengths 1288 INTEGER :: jpimin,jpiwid ! temporary indices 1289 INTEGER :: jpimin_lsm,jpiwid_lsm ! temporary indices 1290 INTEGER :: jpjmin,jpjwid ! temporary indices 1291 INTEGER :: jpjmin_lsm,jpjwid_lsm ! temporary indices 1292 INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices 1293 INTEGER :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices 1294 INTEGER :: itmpi,itmpj,itmpz ! lengths 1295 REAL(wp),DIMENSION(:,:,:), ALLOCATABLE :: ztmp_fly_dta ! local array of values on input grid 1305 1296 !!---------------------------------------------------------------------- 1306 1297 ! … … 1355 1346 1356 1347 1357 itmpi= SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),1)1358 itmpj= SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),2)1348 itmpi=jpi2_lsm-jpi1_lsm+1 1349 itmpj=jpj2_lsm-jpj1_lsm+1 1359 1350 itmpz=kk 1360 1351 ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
Note: See TracChangeset
for help on using the changeset viewer.