Changeset 2051 for branches/DEV_r1784_3DF/NEMO/OPA_SRC/SBC/fldread.F90
- Timestamp:
- 2010-08-13T10:47:35+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r1784_3DF/NEMO/OPA_SRC/SBC/fldread.F90
r1856 r2051 48 48 INTEGER , DIMENSION(2) :: nrec_b ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year) 49 49 INTEGER , DIMENSION(2) :: nrec_a ! after record (1: index, 2: second since Jan. 1st 00h of nit000 year) 50 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,: ):: fnow ! input fields interpolated to now time step50 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:) :: fnow ! input fields interpolated to now time step 51 51 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) :: fdta ! 2 consecutive record of input fields 52 52 CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file acting as a key … … 78 78 INTEGER, DIMENSION(:,:,:), POINTER :: data_jpj ! array of source integers 79 79 REAL(wp), DIMENSION(:,:,:), POINTER :: data_wgt ! array of weights on model grid 80 REAL(wp), DIMENSION(:,: ), POINTER:: fly_dta ! array of values on input grid81 REAL(wp), DIMENSION(:,: ), POINTER:: col2 ! temporary array for reading in columns80 REAL(wp), DIMENSION(:,:,:), POINTER :: fly_dta ! array of values on input grid 81 REAL(wp), DIMENSION(:,:,:), POINTER :: col2 ! temporary array for reading in columns 82 82 END TYPE WGT 83 83 … … 159 159 160 160 ! last record to be read in the current file 161 IF( sd(jf)%nfreqh == -1 ) THEN ; ireclast = 12 161 IF( sd(jf)%nfreqh == -1 ) THEN 162 IF( sd(jf)%cltype == 'monthly' ) THEN ; ireclast = 1 163 ELSE ; ireclast = 12 164 ENDIF 162 165 ELSE 163 166 IF( sd(jf)%cltype == 'monthly' ) THEN ; ireclast = 24 * nmonth_len(nmonth) / sd(jf)%nfreqh … … 204 207 205 208 ! read after data 206 ipk = SIZE( sd(jf)%fdta, 3 )207 209 IF( LEN(TRIM(sd(jf)%wgtname)) > 0 ) THEN 208 210 CALL wgt_list( sd(jf), kw ) 209 DO jk = 1, ipk 210 CALL fld_interp( sd(jf)%num, sd(jf)%clvar, kw, sd(jf)%fdta(:,:,jk,2), sd(jf)%nrec_a(1) ) 211 END DO 211 ipk = SIZE(sd(jf)%fdta,3) 212 CALL fld_interp( sd(jf)%num, sd(jf)%clvar , kw , ipk, sd(jf)%fdta(:,:,:,2) , sd(jf)%nrec_a(1) ) 212 213 ELSE 213 IF( ipk == 1 ) THEN 214 SELECT CASE( SIZE(sd(jf)%fdta,3) ) 215 CASE(1) 214 216 CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fdta(:,:,1,2), sd(jf)%nrec_a(1) ) 215 ELSE217 CASE(jpk) 216 218 CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fdta(:,:,:,2), sd(jf)%nrec_a(1) ) 217 END IF219 END SELECT 218 220 ENDIF 219 221 sd(jf)%rotn(2) = .FALSE. … … 255 257 vtmp(:,:) = 0.0 256 258 ! 257 DO jk = 1, SIZE( sd(kf)%fdta, 3 ) 259 ipk = SIZE( sd(kf)%fdta(:,:,:,nf) ,3 ) 260 DO jk = 1,ipk 258 261 CALL rot_rep( sd(jf)%fdta(:,:,jk,nf),sd(kf)%fdta(:,:,jk,nf),'T', 'en->i', utmp(:,:) ) 259 262 CALL rot_rep( sd(jf)%fdta(:,:,jk,nf),sd(kf)%fdta(:,:,jk,nf),'T', 'en->j', vtmp(:,:) ) … … 333 336 INTEGER :: inrec ! number of record existing for this variable 334 337 INTEGER :: kwgt 335 INTEGER :: jk ! 336 INTEGER :: ipk ! 338 INTEGER :: jk !vertical loop variable 339 INTEGER :: ipk !number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 337 340 CHARACTER(LEN=1000) :: clfmt ! write format 338 341 !!--------------------------------------------------------------------- … … 381 384 & nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /)), & 382 385 & nday - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)), .NOT. llprev ) 383 386 384 387 ! if previous year/month/day file does not exist, we switch to the current year/month/day 385 388 IF( llprev .AND. sdjf%num == 0 ) THEN … … 399 402 400 403 ! read before data into sdjf%fdta(:,:,2) because we will swap data in the following part of fld_read 401 ipk = SIZE( sdjf%fdta, 3 )404 402 405 IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 403 406 CALL wgt_list( sdjf, kwgt ) 404 DO jk = 1, ipk 405 CALL fld_interp( sdjf%num, sdjf%clvar, kwgt, sdjf%fdta(:,:,jk,2), sdjf%nrec_b(1) ) 406 END DO 407 ipk = SIZE(sdjf%fdta,3) 408 CALL fld_interp( sdjf%num, sdjf%clvar, kwgt, ipk, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) ) 407 409 ELSE 408 IF( ipk == 1 ) THEN 409 CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_b(1) ) 410 ELSE 411 CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_b(1) ) 412 ENDIF 410 SELECT CASE ( SIZE(sdjf%fdta,3) ) 411 CASE(1) 412 CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_b(1) ) 413 CASE(jpk) 414 CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_b(1) ) 415 END SELECT 413 416 ENDIF 414 417 sdjf%rotn(2) = .FALSE. … … 421 424 ENDIF 422 425 426 423 427 IF( sdjf%num == 0 ) CALL fld_clopn( sdjf, nyear, nmonth, nday ) ! make sure current year/month/day file is opened 424 428 425 429 sdjf%nswap_sec = nsec_year + nsec1jan000 - 1 ! force read/update the after data in the following part of fld_read 426 430 427 431 END SUBROUTINE fld_init 428 432 … … 458 462 ! forcing record : nmonth 459 463 ! 460 ztmp = REAL( nday, wp ) / REAL( nmonth_len(nmonth), wp ) + 0.5 464 ztmp = 0.e0 465 IF( REAL( nday, wp ) / REAL( nmonth_len(nmonth), wp ) .GT. 0.5 ) ztmp = 1.0 461 466 ELSE 462 467 ztmp = 0.e0 … … 468 473 ENDIF 469 474 470 sdjf%nrec_a(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /) ! define after record number and time 471 irec = irec - 1 ! move back to previous record 472 sdjf%nrec_b(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /) ! define before record number and time 475 IF( sdjf%cltype == 'monthly' ) THEN 476 477 sdjf%nrec_b(:) = (/ 0, nmonth_half(irec - 1 ) + nsec1jan000 /) 478 sdjf%nrec_a(:) = (/ 1, nmonth_half(irec ) + nsec1jan000 /) 479 480 IF( ztmp == 1. ) THEN 481 sdjf%nrec_b(1) = 1 482 sdjf%nrec_a(1) = 2 483 ENDIF 484 485 ELSE 486 487 sdjf%nrec_a(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /) ! define after record number and time 488 irec = irec - 1 ! move back to previous record 489 sdjf%nrec_b(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /) ! define before record number and time 490 491 ENDIF 473 492 ! 474 493 ELSE ! higher frequency mean (in hours) … … 558 577 ELSE 559 578 ! build the new filename if climatological data 560 IF( sdjf%cltype == 'monthly' ) WRITE(sdjf%clname, '(a,"_m" ,i2.2)' ) TRIM( sdjf%clrootname ), kmonth ! add month579 IF( sdjf%cltype == 'monthly' ) WRITE(sdjf%clname, '(a,"_m" ,i2.2)' ) TRIM( sdjf%clrootname ), kmonth ! add month 561 580 ENDIF 562 581 CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 ) … … 707 726 INTEGER :: inum ! temporary logical unit 708 727 INTEGER :: id ! temporary variable id 728 INTEGER :: ipk ! temporary vertical dimension 709 729 CHARACTER (len=5) :: aname 710 730 INTEGER , DIMENSION(3) :: ddims … … 871 891 ! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration. 872 892 ! a more robust solution will be given in next release 873 ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3) ) 874 IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col2(2,ref_wgts(nxt_wgt)%jpjwgt+3) ) 893 ipk = SIZE(sd%fdta,3) 894 ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) ) 895 IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col2(2,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) ) 875 896 876 897 nxt_wgt = nxt_wgt + 1 … … 882 903 END SUBROUTINE fld_weight 883 904 884 SUBROUTINE fld_interp(num, clvar, kw, dta, nrec)905 SUBROUTINE fld_interp(num, clvar, kw, kk, dta, nrec) 885 906 !!--------------------------------------------------------------------- 886 907 !! *** ROUTINE fld_interp *** … … 891 912 !! ** Method : 892 913 !!---------------------------------------------------------------------- 893 INTEGER, INTENT(in) :: num ! stream number 894 CHARACTER(LEN=*), INTENT(in) :: clvar ! variable name 895 INTEGER, INTENT(in) :: kw ! weights number 896 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: dta ! output field on model grid 897 INTEGER, INTENT(in) :: nrec ! record number to read (ie time slice) 914 INTEGER, INTENT(in) :: num ! stream number 915 CHARACTER(LEN=*), INTENT(in) :: clvar ! variable name 916 INTEGER, INTENT(in) :: kw ! weights number 917 INTEGER, INTENT(in) :: kk ! vertical dimension of kk 918 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kk) :: dta ! output field on model grid 919 INTEGER, INTENT(in) :: nrec ! record number to read (ie time slice) 898 920 !! 899 INTEGER, DIMENSION( 2):: rec1,recn ! temporary arrays for start and length900 INTEGER :: jk, jn, jm ! loop counters901 INTEGER :: ni, nj ! lengths902 INTEGER :: jpimin,jpiwid ! temporary indices903 INTEGER :: jpjmin,jpjwid ! temporary indices904 INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices921 INTEGER, DIMENSION(3) :: rec1,recn ! temporary arrays for start and length 922 INTEGER :: jk, jn, jm ! loop counters 923 INTEGER :: ni, nj ! lengths 924 INTEGER :: jpimin,jpiwid ! temporary indices 925 INTEGER :: jpjmin,jpjwid ! temporary indices 926 INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices 905 927 !!---------------------------------------------------------------------- 906 928 ! … … 920 942 rec1(1) = MAX( jpimin-1, 1 ) 921 943 rec1(2) = MAX( jpjmin-1, 1 ) 944 rec1(3) = 1 922 945 recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 ) 923 946 recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 ) 947 recn(3) = kk 924 948 925 949 !! where we need to read it to … … 929 953 jpj2 = jpj1 + recn(2) - 1 930 954 931 ref_wgts(kw)%fly_dta(:,:) = 0.0 932 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2), nrec, rec1, recn) 955 ref_wgts(kw)%fly_dta(:,:,:) = 0.0 956 SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) ) 957 CASE(1) 958 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn) 959 CASE(jpk) 960 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn) 961 END SELECT 933 962 934 963 !! first four weights common to both bilinear and bicubic 935 964 !! note that we have to offset by 1 into fly_dta array because of halo 936 dta(:,: ) = 0.0965 dta(:,:,:) = 0.0 937 966 DO jk = 1,4 938 DO jn = 1, jpj939 DO jm = 1, jpi967 DO jn = 1, nlcj 968 DO jm = 1,nlci 940 969 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 941 970 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 942 dta(jm,jn ) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1)971 dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,jk) 943 972 END DO 944 973 END DO … … 949 978 !! fix up halo points that we couldnt read from file 950 979 IF( jpi1 == 2 ) THEN 951 ref_wgts(kw)%fly_dta(jpi1-1,: ) = ref_wgts(kw)%fly_dta(jpi1,:)980 ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:) 952 981 ENDIF 953 982 IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 954 ref_wgts(kw)%fly_dta(jpi2+1,: ) = ref_wgts(kw)%fly_dta(jpi2,:)983 ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:) 955 984 ENDIF 956 985 IF( jpj1 == 2 ) THEN 957 ref_wgts(kw)%fly_dta(:,jpj1-1 ) = ref_wgts(kw)%fly_dta(:,jpj1)986 ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:) 958 987 ENDIF 959 988 IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN 960 ref_wgts(kw)%fly_dta(:,jpj2+1 ) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2) - ref_wgts(kw)%fly_dta(:,jpj2-1)989 ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:) 961 990 ENDIF 962 991 … … 971 1000 IF( jpi1 == 2 ) THEN 972 1001 rec1(1) = ref_wgts(kw)%ddims(1) - 1 973 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn) 974 ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2) = ref_wgts(kw)%col2(ref_wgts(kw)%offset+1,jpj1:jpj2) 1002 SELECT CASE( SIZE( ref_wgts(kw)%col2(:,jpj1:jpj2,:),3) ) 1003 CASE(1) 1004 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,1), nrec, rec1, recn) 1005 CASE(jpk) 1006 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,:), nrec, rec1, recn) 1007 END SELECT 1008 ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col2(ref_wgts(kw)%offset+1,jpj1:jpj2,:) 975 1009 ENDIF 976 1010 IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 977 1011 rec1(1) = 1 978 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn) 979 ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2) = ref_wgts(kw)%col2(2-ref_wgts(kw)%offset,jpj1:jpj2) 1012 SELECT CASE( SIZE( ref_wgts(kw)%col2(:,jpj1:jpj2,:),3) ) 1013 CASE(1) 1014 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,1), nrec, rec1, recn) 1015 CASE(jpk) 1016 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,:), nrec, rec1, recn) 1017 END SELECT 1018 ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col2(2-ref_wgts(kw)%offset,jpj1:jpj2,:) 980 1019 ENDIF 981 1020 ENDIF … … 983 1022 ! gradient in the i direction 984 1023 DO jk = 1,4 985 DO jn = 1, jpj986 DO jm = 1, jpi1024 DO jn = 1, nlcj 1025 DO jm = 1,nlci 987 1026 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 988 1027 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 989 dta(jm,jn ) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * &990 (ref_wgts(kw)%fly_dta(ni+2,nj+1 ) - ref_wgts(kw)%fly_dta(ni,nj+1))1028 dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * & 1029 (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:)) 991 1030 END DO 992 1031 END DO … … 995 1034 ! gradient in the j direction 996 1035 DO jk = 1,4 997 DO jn = 1, jpj998 DO jm = 1, jpi1036 DO jn = 1, nlcj 1037 DO jm = 1,nlci 999 1038 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 1000 1039 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 1001 dta(jm,jn ) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * &1002 (ref_wgts(kw)%fly_dta(ni+1,nj+2 ) - ref_wgts(kw)%fly_dta(ni+1,nj))1040 dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * & 1041 (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:)) 1003 1042 END DO 1004 1043 END DO … … 1011 1050 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 1012 1051 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 1013 dta(jm,jn ) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &1014 (ref_wgts(kw)%fly_dta(ni+2,nj+2 ) - ref_wgts(kw)%fly_dta(ni ,nj+2)) - &1015 (ref_wgts(kw)%fly_dta(ni+2,nj ) - ref_wgts(kw)%fly_dta(ni ,nj)))1052 dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( & 1053 (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni ,nj+2,:)) - & 1054 (ref_wgts(kw)%fly_dta(ni+2,nj ,:) - ref_wgts(kw)%fly_dta(ni ,nj ,:))) 1016 1055 END DO 1017 1056 END DO
Note: See TracChangeset
for help on using the changeset viewer.
