From 7e93992066eaece4f004cc5392a8ebe0f8c16717 Mon Sep 17 00:00:00 2001 From: pdicerbo Date: Thu, 27 Oct 2016 11:00:47 +0200 Subject: [PATCH] def_cov improvement --- CMakeLists.txt | 2 +- Makefile | 2 + cns_str.f90 | 7 +- def_cov.f90 | 238 +++++++++++++++++++++++++++-------------------- filename_mod.f90 | 4 + get_obs_arg.f90 | 10 +- get_obs_chl.f90 | 4 +- grd_str.f90 | 16 ++-- rcfl_y_init.f90 | 2 +- rdeofs.f90 | 4 +- rdgrds.f90 | 18 ++-- ver_hor.f90 | 8 +- ver_hor_ad.f90 | 8 +- wrt_dia.f90 | 2 +- 14 files changed, 185 insertions(+), 140 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e6710dc..cfac64b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -24,7 +24,7 @@ set(OBJS routines.f def_nml.f90 def_grd.f90 sav_itr.f90 ini_itr.f90 ini_bmd.f90 netcdf_err.f90 get_obs.f90get_obs_chl.f90 int_par.f90 obs_vec.f90 def_cov.f90 ini_cfn.f90 ini_nrm.f90 min_cfn.f90 costf.f90 cnv_ctv.f90 ver_hor.f90 rcfl_x.f90 rcfl_y.f90 veof.f90 - obsop.f90 obs_chl.f90 resid.f90 res_inc.f90 + obsop.f90 obs_chl.f90 resid.f90 res_inc.f90 rdrcorr.f90 mean_rdr.f90 obs_trd_ad.f90 obs_chl_ad.f90 veof_ad.f90 ver_hor_ad.f90 rcfl_x_ad.f90 rcfl_y_ad.f90 cnv_ctv_ad.f90 cnv_inn.f90 wrt_dia.f90 filename_mod.f90 oceanvar.f90 main.f90) # SOURCES --END-- diff --git a/Makefile b/Makefile index 0b28977..93df097 100644 --- a/Makefile +++ b/Makefile @@ -106,6 +106,8 @@ OBJS = \ ini_itr.o\ rdgrds.o\ rdeofs.o\ + rdrcorr.o\ + mean_rdr.o\ netcdf_err.o\ get_obs.o\ get_obs_arg.o\ diff --git a/cns_str.f90 b/cns_str.f90 index 9222961..041ee2d 100644 --- a/cns_str.f90 +++ b/cns_str.f90 @@ -39,6 +39,9 @@ MODULE cns_str INTEGER(i4) :: ntr ! No. of iterations (half of) REAL(r8) :: dx ! Grid resolution (m) REAL(r8) :: L ! Correlation radius +!laura + REAL(r8),POINTER :: Lxyz(:,:,:)!Correlation radius from file in km +!laura REAL(r8) :: E ! Norm REAL(r8) :: alp ! Filter weight INTEGER(i4) :: ntb ! Number of points in the table @@ -46,7 +49,9 @@ MODULE cns_str REAL(r8) :: dsmx ! Maximum distance REAL(r8) :: dsl ! Table increment REAL(r8), POINTER :: al(:) ! Filter weights in the table - REAL(r8), POINTER :: sc(:) ! Filter scaling factors in the table +!laura + REAL(r8), POINTER :: sc(:,:) ! Filter scaling factors in the table +!laura REAL(r8) :: scl ! Scaling factor REAL(r8) :: efc ! Scaling factor for extended points INTEGER(i4) :: kstp ! Step for extended points diff --git a/def_cov.f90 b/def_cov.f90 index 28e7af8..69e34b7 100644 --- a/def_cov.f90 +++ b/def_cov.f90 @@ -37,12 +37,17 @@ subroutine def_cov implicit none - INTEGER(i4) :: k, nspl, i, j, kk - REAL(r8) :: E, dst + INTEGER(i4) :: k, nspl, i, j, kk,l,m + REAL(r8) :: dst,E,Lmean, mean_rad REAL(r8) , ALLOCATABLE :: sfct(:), al(:), bt(:) INTEGER(i4) , ALLOCATABLE :: jnxx(:) INTEGER nthreads, threadid integer :: OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM + + !laura + call rdrcorr !laura Allocation rcf%Lxyz + !laura + nthreads = 1 threadid = 0 @@ -70,9 +75,9 @@ subroutine def_cov ! KB grid problem (chl assimilation) rcf%ntb = 1000 ! - ALLOCATE ( rcf%al(rcf%ntb)) ; rcf%al = huge(rcf%al(1)) - ALLOCATE ( rcf%sc(rcf%ntb)) ; rcf%sc = huge(rcf%sc(1)) + !laura, sc added vertical level + ALLOCATE ( rcf%sc(grd%km,rcf%ntb)) ; rcf%sc = huge(rcf%sc(1,1)) rcf%dsmn = 1.e20 rcf%dsmx = -1.e20 @@ -86,118 +91,144 @@ subroutine def_cov rcf%dsmx = rcf%dsmx + max(1.d0,(rcf%dsmx-rcf%dsmn)/(rcf%ntb-2.)) rcf%dsl = (rcf%dsmx-rcf%dsmn) / (rcf%ntb-1.) - - do k=1,rcf%ntb - dst = rcf%dsmn + (k-1.) * rcf%dsl - E = (2. * rcf%ntr) * dst**2 / (4. * rcf%L**2) - rcf%al(k) = 1. + E - sqrt(E*(E+2.)) - rcf%alp = rcf%al(k) - sfct(:) = 0. - al(:) = rcf%al(k) - bt(:) = rcf%al(k) - do j=1,nspl - jnxx(j) = j - enddo - sfct(nspl/2+1) = 1. - call rcfl_y_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nspl) - call rcfl_y_ad_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nspl) - rcf%sc(k) = sfct(nspl/2+1) - enddo + +!Laura aggiunto ciclo sui livelli + do k=1,grd%km + Lmean=mean_rad(k,rcf%Lxyz(:,:,k)) !for each level, mean of Lxyz + do l=1,rcf%ntb + dst = rcf%dsmn + (l-1.) * rcf%dsl + !E = (2. * rcf%ntr) * dst**2 / (4. * rcf%L**2) + E = (2. * rcf%ntr) * dst**2 / (4. * Lmean**2) !mean of radius + rcf%al(l) = 1. + E - sqrt(E*(E+2.)) + rcf%alp = rcf%al(l) + sfct(:) = 0. + al(:) = rcf%al(l) + bt(:) = rcf%al(l) + do m=1,nspl + jnxx(m) = m + enddo + sfct(nspl/2+1) = 1. + call rcfl_y_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nspl) + call rcfl_y_ad_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nspl) + rcf%sc(k,l) = sfct(nspl/2+1) + enddo + enddo !k vertical loop DEALLOCATE ( sfct, jnxx, al, bt ) - do j=1,grd%jm - do i=1,grd%im - dst = ( grd%dx(i,j) - rcf%dsmn )/rcf%dsl - k = int(dst) + 1 - dst = dst - real(k-1) - grd%scx(i,j) = sqrt( 1./ (rcf%sc(k)*(1.-dst) + rcf%sc(k+1)*dst) ) - dst = ( grd%dy(i,j) - rcf%dsmn )/rcf%dsl - k = int(dst) + 1 - dst = dst - real(k-1) - grd%scy(i,j) = sqrt( 1./ (rcf%sc(k)*(1.-dst) + rcf%sc(k+1)*dst) ) - enddo - enddo - - do j=1,grd%jm - do i=2,grd%im - dst = (grd%dx(i-1,j) + grd%dx(i,j)) * 0.5 - E = (2. * rcf%ntr) * dst**2 / (4. * rcf%L**2) - grd%alx(i,j) = 1. + E - sqrt(E*(E+2.)) - enddo - do i=1,grd%im-1 - dst = (grd%dx(i,j) + grd%dx(i+1,j)) * 0.5 - E = (2. * rcf%ntr) * dst**2 / (4. * rcf%L**2) - grd%btx(i,j) = 1. + E - sqrt(E*(E+2.)) - enddo - enddo - - do j=2,grd%jm - do i=1,grd%im - dst = (grd%dy(i,j-1) + grd%dy(i,j)) * 0.5 - E = (2. * rcf%ntr) * dst**2 / (4. * rcf%L**2) - grd%aly(i,j) = 1. + E - sqrt(E*(E+2.)) - enddo - enddo - do j=1,grd%jm-1 - do i=1,grd%im - dst = (grd%dy(i,j) + grd%dy(i,j+1)) * 0.5 - E = (2. * rcf%ntr) * dst**2 / (4. * rcf%L**2) - grd%bty(i,j) = 1. + E - sqrt(E*(E+2.)) - enddo - enddo - - grd%alx( 1,:) = grd%alx( 2,:) - grd%btx(grd%im,:) = grd%btx(grd%im-1,:) - grd%aly(:, 1) = grd%aly(:, 2) - grd%bty(:,grd%jm) = grd%bty(:,grd%jm-1) - - !--- - ! Define extended grids - - grd%istp = int( rcf%L * rcf%efc / grd%dx(:,:) )+1 - grd%jstp = int( rcf%L * rcf%efc / grd%dy(:,:) )+1 +! grd%scx=1.e20 !add undef values for scx and scy +! grd%scy=1.e20 + do k=1,grd%km + do j=1,grd%jm + do i=1,grd%im + dst = ( grd%dx(i,j) - rcf%dsmn )/rcf%dsl + l = int(dst) + 1 + dst = dst - real(l-1) + grd%scx(i,j,k) = sqrt( 1./ (rcf%sc(k,l)*(1.-dst) + rcf%sc(k,l+1)*dst) ) + dst = ( grd%dy(i,j) - rcf%dsmn )/rcf%dsl + l = int(dst) + 1 + dst = dst - real(l-1) + grd%scy(i,j,k) = sqrt( 1./ (rcf%sc(k,l)*(1.-dst) + rcf%sc(k,l+1)*dst) ) + enddo + enddo + enddo !k vertical loop + + +! modificato fin qui + + + do k=1,grd%km !laura add vertical component + do j=1,grd%jm + do i=2,grd%im + dst = (grd%dx(i-1,j) + grd%dx(i,j)) * 0.5 + E = (2. * rcf%ntr) * dst**2 / (4. * rcf%Lxyz(i,j,k)**2) !laura + grd%alx(i,j,k) = 1. + E - sqrt(E*(E+2.)) + enddo + do i=1,grd%im-1 + dst = (grd%dx(i,j) + grd%dx(i+1,j)) * 0.5 + E = (2. * rcf%ntr) * dst**2 / (4. * rcf%Lxyz(i,j,k)**2) !laura + grd%btx(i,j,k) = 1. + E - sqrt(E*(E+2.)) + enddo + enddo + + do j=2,grd%jm + do i=1,grd%im + dst = (grd%dy(i,j-1) + grd%dy(i,j)) * 0.5 + E = (2. * rcf%ntr) * dst**2 / (4. * rcf%Lxyz(i,j,k)**2) !laura + grd%aly(i,j,k) = 1. + E - sqrt(E*(E+2.)) + enddo + enddo + do j=1,grd%jm-1 + do i=1,grd%im + dst = (grd%dy(i,j) + grd%dy(i,j+1)) * 0.5 + E = (2. * rcf%ntr) * dst**2 / (4. * rcf%Lxyz(i,j,k)**2) !laura + grd%bty(i,j,k) = 1. + E - sqrt(E*(E+2.)) + enddo + enddo + + grd%alx( 1,:,k) = grd%alx( 2,:,k) + grd%btx(grd%im,:,k) = grd%btx(grd%im-1,:,k) + grd%aly(:, 1,k) = grd%aly(:, 2,k) + grd%bty(:,grd%jm,k) = grd%bty(:,grd%jm-1,k) + + !--- + ! Define extended grids + + !grd%istp = int( rcf%L * rcf%efc / grd%dx(:,:) )+1 + !grd%jstp = int( rcf%L * rcf%efc / grd%dy(:,:) )+1 + grd%istp(:,:,k) = int( rcf%Lxyz(:,:,k) * rcf%efc / grd%dx(:,:) )+1 !laura + grd%jstp(:,:,k) = int( rcf%Lxyz(:,:,k) * rcf%efc / grd%dy(:,:) )+1 + enddo ! k vertical component + + grd%imax = 0 grd%jmax = 0 - do k = 1, grd%km grd%imx(k) = 0 do j = 1, grd%jm - kk = grd%istp(1,j) + !kk = grd%istp(1,j) + kk = grd%istp(1,j,k) !laura aggiunta la verticale if( grd%msr(1,j,k).eq.1. ) kk = kk + 1 grd%inx(1,j,k) = kk do i = 2, grd%im if( grd%msr(i,j,k).eq.0. .and. grd%msr(i-1,j,k).eq.1. ) then - kk = kk + grd%istp(i,j) + !kk = kk + grd%istp(i,j) + kk = kk + grd%istp(i,j,k) !laura else if( grd%msr(i,j,k).eq.1. .and. grd%msr(i-1,j,k).eq.0. ) then - kk = kk + grd%istp(i,j) + 1 + !kk = kk + grd%istp(i,j) + 1 + kk = kk + grd%istp(i,j,k) + 1 !laura else if( grd%msr(i,j,k).eq.1. ) then kk = kk + 1 endif grd%inx(i,j,k) = kk enddo - grd%imx(k) = max( grd%imx(k), kk+grd%istp(grd%im,j)) + !grd%imx(k) = max( grd%imx(k), kk+grd%istp(grd%im,j)) + grd%imx(k) = max( grd%imx(k), kk+grd%istp(grd%im,j,k)) !laura enddo grd%imax = max( grd%imax, grd%imx(k)) grd%jmx(k) = 0 do i = 1, grd%im - kk = grd%jstp(i,1) + !kk = grd%jstp(i,1) + kk = grd%jstp(i,1,k) !laura if( grd%msr(i,1,k).eq.1. ) kk = kk + 1 grd%jnx(i,1,k) = kk do j = 2, grd%jm if( grd%msr(i,j,k).eq.0. .and. grd%msr(i,j-1,k).eq.1. ) then - kk = kk + grd%jstp(i,j) + !kk = kk + grd%jstp(i,j) + kk = kk + grd%jstp(i,j,k) !laura else if( grd%msr(i,j,k).eq.1. .and. grd%msr(i,j-1,k).eq.0. ) then - kk = kk + grd%jstp(i,j) + 1 + !kk = kk + grd%jstp(i,j) + 1 + kk = kk + grd%jstp(i,j,k) + 1 !laura else if( grd%msr(i,j,k).eq.1. ) then kk = kk + 1 endif grd%jnx(i,j,k) = kk enddo - grd%jmx(k) = max( grd%jmx(k), kk+grd%jstp(i,grd%jm)) + !grd%jmx(k) = max( grd%jmx(k), kk+grd%jstp(i,grd%jm)) + grd%jmx(k) = max( grd%jmx(k), kk+grd%jstp(i,grd%jm,k)) !laura enddo grd%jmax = max( grd%jmax, grd%jmx(k)) @@ -213,48 +244,49 @@ subroutine def_cov do k = 1, grd%km do j = 1, grd%jm - kk = grd%istp(1,j) + !kk = grd%istp(1,j) + kk = grd%istp(1,j,k) !laura if( grd%msr(1,j,k).eq.1. ) then kk = kk + 1 - grd%aex(j,1:kk,k) = grd%alx(1,j) - grd%bex(j,1:kk,k) = grd%btx(1,j) + grd%aex(j,1:kk,k) = grd%alx(1,j,k) + grd%bex(j,1:kk,k) = grd%btx(1,j,k) endif do i = 2, grd%im if( grd%msr(i,j,k).eq.0. .and. grd%msr(i-1,j,k).eq.1. ) then - grd%aex(j,kk+1:kk+grd%istp(i,j),k) = grd%alx(i,j) - grd%bex(j,kk+1:kk+grd%istp(i,j),k) = grd%btx(i,j) - kk = kk + grd%istp(i,j) + grd%aex(j,kk+1:kk+grd%istp(i,j,k),k) = grd%alx(i,j,k) + grd%bex(j,kk+1:kk+grd%istp(i,j,k),k) = grd%btx(i,j,k) + kk = kk + grd%istp(i,j,k) else if( grd%msr(i,j,k).eq.1. .and. grd%msr(i-1,j,k).eq.0. ) then - grd%aex(j,kk+1:kk+grd%istp(i,j)+1,k) = grd%alx(i,j) - grd%bex(j,kk+1:kk+grd%istp(i,j)+1,k) = grd%btx(i,j) - kk = kk + grd%istp(i,j) + 1 + grd%aex(j,kk+1:kk+grd%istp(i,j,k)+1,k) = grd%alx(i,j,k) + grd%bex(j,kk+1:kk+grd%istp(i,j,k)+1,k) = grd%btx(i,j,k) + kk = kk + grd%istp(i,j,k) + 1 else if( grd%msr(i,j,k).eq.1. ) then - grd%aex(j,kk+1,k) = grd%alx(i,j) - grd%bex(j,kk+1,k) = grd%btx(i,j) + grd%aex(j,kk+1,k) = grd%alx(i,j,k) + grd%bex(j,kk+1,k) = grd%btx(i,j,k) kk = kk + 1 endif enddo enddo do i = 1, grd%im - kk = grd%jstp(i,1) + kk = grd%jstp(i,1,k) if( grd%msr(i,1,k).eq.1. ) then kk = kk + 1 - grd%aey(i,1:kk,k) = grd%aly(i,1) - grd%bey(i,1:kk,k) = grd%bty(i,1) + grd%aey(i,1:kk,k) = grd%aly(i,1,k) + grd%bey(i,1:kk,k) = grd%bty(i,1,k) endif do j = 2, grd%jm if( grd%msr(i,j,k).eq.0. .and. grd%msr(i,j-1,k).eq.1. ) then - grd%aey(i,kk+1:kk+grd%jstp(i,j),k) = grd%aly(i,j) - grd%bey(i,kk+1:kk+grd%jstp(i,j),k) = grd%bty(i,j) - kk = kk + grd%jstp(i,j) + grd%aey(i,kk+1:kk+grd%jstp(i,j,k),k) = grd%aly(i,j,k) + grd%bey(i,kk+1:kk+grd%jstp(i,j,k),k) = grd%bty(i,j,k) + kk = kk + grd%jstp(i,j,k) else if( grd%msr(i,j,k).eq.1. .and. grd%msr(i,j-1,k).eq.0. ) then - grd%aey(i,kk+1:kk+grd%jstp(i,j)+1,k) = grd%aly(i,j) - grd%bey(i,kk+1:kk+grd%jstp(i,j)+1,k) = grd%bty(i,j) - kk = kk + grd%jstp(i,j) + 1 + grd%aey(i,kk+1:kk+grd%jstp(i,j,k)+1,k) = grd%aly(i,j,k) + grd%bey(i,kk+1:kk+grd%jstp(i,j,k)+1,k) = grd%bty(i,j,k) + kk = kk + grd%jstp(i,j,k) + 1 else if( grd%msr(i,j,k).eq.1. ) then - grd%aey(i,kk+1,k) = grd%aly(i,j) - grd%bey(i,kk+1,k) = grd%bty(i,j) + grd%aey(i,kk+1,k) = grd%aly(i,j,k) + grd%bey(i,kk+1,k) = grd%bty(i,j,k) kk = kk + 1 endif enddo diff --git a/filename_mod.f90 b/filename_mod.f90 index e9da2e3..8040003 100644 --- a/filename_mod.f90 +++ b/filename_mod.f90 @@ -8,6 +8,8 @@ MODULE FILENAMES character (LEN=1024) :: EIV_FILE != 'eiv.nc' character (LEN=1024) :: OBS_FILE != 'obs_1.dat' character (LEN=1024) :: GRID_FILE != 'grid1.nc' +!laura +character (LEN=1024) :: RCORR_FILE != 'chl_rad_corr.nc' @@ -24,6 +26,8 @@ SUBROUTINE SETFILENAMES EIV_FILE = 'eiv.nc' OBS_FILE = 'obs_1.dat' ! 'obs_'//fgrd//'.dat' GRID_FILE = 'grid1.nc'! 'grid'//cgrd//'.nc' +!laura + RCORR_FILE = 'chl_rad_corr.nc' END SUBROUTINE SETFILENAMES diff --git a/get_obs_arg.f90 b/get_obs_arg.f90 index 529e9ee..1c1e498 100644 --- a/get_obs_arg.f90 +++ b/get_obs_arg.f90 @@ -42,7 +42,7 @@ subroutine get_obs_arg arg%nc = 0 - open(511,file='arg_mis.dat',form='formatted') + open(511,file='arg_datnew.dat',form='formatted') ! --- ! Allocate memory for observations @@ -67,7 +67,7 @@ subroutine get_obs_arg arg%bia(:) = 0.0 do k=1,arg%no - read (511,'(I5,I5,F10.3,F10.3,F10.3,F10.3,F10.3,F10.3,I8)') & + read (511,'(I5,I5,F12.5,F12.5,F12.5,F12.5,F12.5,F12.5,I8)') & arg%flc(k), arg%par(k), & arg%lon(k), arg%lat(k), & arg%dpt(k), arg%tim(k), & @@ -75,9 +75,10 @@ subroutine get_obs_arg end do close (511) + ! DECOMMENT FOLLOWING TWO LINES TO MAKE FILTER TEST - arg%res(:) = 1 - arg%err(:) = 1.d-2 + !arg%res(:) = 1 + !arg%err(:) = 1.d-2 ! --- @@ -125,7 +126,6 @@ subroutine get_obs_arg arg%pq8(k) = 0. endif enddo - arg%flc(:) = arg%flg(:) end subroutine get_obs_arg diff --git a/get_obs_chl.f90 b/get_obs_chl.f90 index a60c070..10f329f 100644 --- a/get_obs_chl.f90 +++ b/get_obs_chl.f90 @@ -122,8 +122,8 @@ subroutine get_obs_chl enddo ! DECOMMENT FOLLOWING TWO LINES TO MAKE FILTER TEST - chl%res(:) = 0. - chl%err(:) = 1. + ! chl%res(:) = 0. + ! chl%err(:) = 1. DEALLOCATE( chl_mis ) DEALLOCATE( chl_err ) diff --git a/grd_str.f90 b/grd_str.f90 index dc2cf05..1adb1a7 100644 --- a/grd_str.f90 +++ b/grd_str.f90 @@ -59,19 +59,19 @@ MODULE grd_str REAL(r8), POINTER :: dy(:,:) ! dy - REAL(r8), POINTER :: alx(:,:) ! Coefficient for the positive direction of the recursive filter - REAL(r8), POINTER :: aly(:,:) ! Coefficient for the positive direction of the recursive filter - REAL(r8), POINTER :: btx(:,:) ! Coefficient for the negative direction of the recursive filter - REAL(r8), POINTER :: bty(:,:) ! Coefficient for the negative direction of the recursive filter - REAL(r8), POINTER :: scx(:,:) ! Scaling factor for x direction - REAL(r8), POINTER :: scy(:,:) ! Scaling factor for y direction + REAL(r8), POINTER :: alx(:,:,:) ! Coefficient for the positive direction of the recursive filter + REAL(r8), POINTER :: aly(:,:,:) ! Coefficient for the positive direction of the recursive filter + REAL(r8), POINTER :: btx(:,:,:) ! Coefficient for the negative direction of the recursive filter + REAL(r8), POINTER :: bty(:,:,:) ! Coefficient for the negative direction of the recursive filter + REAL(r8), POINTER :: scx(:,:,:) ! Scaling factor for x direction !laura + REAL(r8), POINTER :: scy(:,:,:) ! Scaling factor for y direction !laura REAL(r8), POINTER :: msr(:,:,:) ! Sea-land mask used in the recursive filter INTEGER(i4) :: imax ! Maximum number of extended points INTEGER(i4) :: jmax ! Maximum number of extended points INTEGER(i4), POINTER :: imx(:) ! Max. no. of extended pnts at each level INTEGER(i4), POINTER :: jmx(:) ! Max. no. of extended pnts at each level - INTEGER(i4), POINTER :: istp(:,:) ! Extended points - INTEGER(i4), POINTER :: jstp(:,:) ! Extended points + INTEGER(i4), POINTER :: istp(:,:,:) ! Extended points + INTEGER(i4), POINTER :: jstp(:,:,:) ! Extended points INTEGER(i4), POINTER :: inx(:,:,:) ! Pointer for extended grid INTEGER(i4), POINTER :: jnx(:,:,:) ! Pointer for extended grid REAL(r8), POINTER :: fct(:,:,:) ! Normalisation factor diff --git a/rcfl_y_init.f90 b/rcfl_y_init.f90 index 2017b86..a24a803 100644 --- a/rcfl_y_init.f90 +++ b/rcfl_y_init.f90 @@ -62,7 +62,7 @@ subroutine rcfl_y_init( im, jm, km, jmax, al, bt, fld, jnx, jmx) bta(:,:) = bt(:,:,k) - do ktr = 1,rcf%ntr + do ktr = 1,rcf%ntr !numero di cicli del filtro ! positive direction if( ktr.eq.1 )then diff --git a/rdeofs.f90 b/rdeofs.f90 index ce02b7c..f2abbeb 100644 --- a/rdeofs.f90 +++ b/rdeofs.f90 @@ -100,8 +100,8 @@ subroutine rdeofs if (stat /= nf90_noerr) call netcdf_err(stat) ! DECOMMENT FOLLOWING TWO LINES TO MAKE FILTER TEST - ros%evc(:,:,:) = 1. - ros%eva(:,:) = 1. + ! ros%evc(:,:,:) = 1. + ! ros%eva(:,:) = 1. stat = nf90_close(ncid) diff --git a/rdgrds.f90 b/rdgrds.f90 index 3e85562..056bdcc 100644 --- a/rdgrds.f90 +++ b/rdgrds.f90 @@ -72,17 +72,19 @@ subroutine rdgrd ALLOCATE ( grd%dx(grd%im,grd%jm)) ; grd%dx = huge(grd%dx(1,1)) ALLOCATE ( grd%dy(grd%im,grd%jm)) ; grd%dy = huge(grd%dy(1,1)) - ALLOCATE ( grd%alx(grd%im,grd%jm) ) ; grd%alx = huge(grd%alx(1,1)) - ALLOCATE ( grd%aly(grd%im,grd%jm) ) ; grd%aly = huge(grd%aly(1,1)) - ALLOCATE ( grd%btx(grd%im,grd%jm) ) ; grd%btx = huge(grd%btx(1,1)) - ALLOCATE ( grd%bty(grd%im,grd%jm) ) ; grd%bty = huge(grd%bty(1,1)) - ALLOCATE ( grd%scx(grd%im,grd%jm) ) ; grd%scx = huge(grd%scx(1,1)) - ALLOCATE ( grd%scy(grd%im,grd%jm) ) ; grd%scy = huge(grd%scy(1,1)) + ALLOCATE ( grd%alx(grd%im,grd%jm,grd%km) ) ; grd%alx = huge(grd%alx(1,1,1)) !laura + ALLOCATE ( grd%aly(grd%im,grd%jm,grd%km) ) ; grd%aly = huge(grd%aly(1,1,1)) !laura + ALLOCATE ( grd%btx(grd%im,grd%jm,grd%km) ) ; grd%btx = huge(grd%btx(1,1,1)) !laura + ALLOCATE ( grd%bty(grd%im,grd%jm,grd%km) ) ; grd%bty = huge(grd%bty(1,1,1)) !laura + ALLOCATE ( grd%scx(grd%im,grd%jm,grd%km) ) ; grd%scx = huge(grd%scx(1,1,1))!laura + ALLOCATE ( grd%scy(grd%im,grd%jm,grd%km) ) ; grd%scy = huge(grd%scy(1,1,1)) ALLOCATE ( grd%msr(grd%im,grd%jm,grd%km) ) ; grd%msr = huge(grd%msr(1,1,1)) ALLOCATE ( grd%imx(grd%km)) ; grd%imx = huge(grd%imx(1)) ALLOCATE ( grd%jmx(grd%km)) ; grd%jmx = huge(grd%jmx(1)) - ALLOCATE ( grd%istp(grd%im,grd%jm)) ; grd%istp = huge(grd%istp(1,1)) - ALLOCATE ( grd%jstp(grd%im,grd%jm)) ; grd%jstp = huge(grd%jstp(1,1)) +!laura--- + ALLOCATE ( grd%istp(grd%im,grd%jm,grd%km)) ; grd%istp = huge(grd%istp(1,1,1)) + ALLOCATE ( grd%jstp(grd%im,grd%jm,grd%km)) ; grd%jstp = huge(grd%jstp(1,1,1)) +!laura --- ALLOCATE ( grd%inx(grd%im,grd%jm,grd%km)) ; grd%inx = huge(grd%inx(1,1,1)) ALLOCATE ( grd%jnx(grd%im,grd%jm,grd%km)) ; grd%jnx = huge(grd%jnx(1,1,1)) ALLOCATE ( grd%fct(grd%im,grd%jm,grd%km) ) ; grd%fct = huge(grd%fct(1,1,1)) diff --git a/ver_hor.f90 b/ver_hor.f90 index ef64f2e..484f59d 100644 --- a/ver_hor.f90 +++ b/ver_hor.f90 @@ -112,7 +112,7 @@ subroutine ver_hor !$OMP PRIVATE(k) !$OMP DO do k=1,grd%km - grd%chl(:,:,k,l) = grd%chl(:,:,k,l) * grd%scx(:,:) + grd%chl(:,:,k,l) = grd%chl(:,:,k,l) * grd%scx(:,:,k) !laura enddo !$OMP END DO !$OMP END PARALLEL @@ -129,7 +129,7 @@ subroutine ver_hor !$OMP PRIVATE(k) !$OMP DO do k=1,grd%km - grd%chl(:,:,k,l) = grd%chl(:,:,k,l) * grd%scy(:,:) + grd%chl(:,:,k,l) = grd%chl(:,:,k,l) * grd%scy(:,:,k) !laura enddo !$OMP END DO !$OMP END PARALLEL @@ -146,7 +146,7 @@ subroutine ver_hor !$OMP PRIVATE(k) !$OMP DO do k=1,grd%km - grd%chl_ad(:,:,k,l) = grd%chl_ad(:,:,k,l) * grd%scy(:,:) + grd%chl_ad(:,:,k,l) = grd%chl_ad(:,:,k,l) * grd%scy(:,:,k) !laura enddo !$OMP END DO !$OMP END PARALLEL @@ -163,7 +163,7 @@ subroutine ver_hor !$OMP PRIVATE(k) !$OMP DO do k=1,grd%km - grd%chl_ad(:,:,k,l) = grd%chl_ad(:,:,k,l) * grd%scx(:,:) + grd%chl_ad(:,:,k,l) = grd%chl_ad(:,:,k,l) * grd%scx(:,:,k) !laura enddo !$OMP END DO !$OMP END PARALLEL diff --git a/ver_hor_ad.f90 b/ver_hor_ad.f90 index d16a8c5..8f6f551 100644 --- a/ver_hor_ad.f90 +++ b/ver_hor_ad.f90 @@ -103,7 +103,7 @@ subroutine ver_hor_ad !$OMP PRIVATE(k) !$OMP DO do k=1,grd%km - grd%chl(:,:,k,l) = grd%chl(:,:,k,l) * grd%scx(:,:) + grd%chl(:,:,k,l) = grd%chl(:,:,k,l) * grd%scx(:,:,k) !laura enddo !$OMP END DO !$OMP END PARALLEL @@ -120,7 +120,7 @@ subroutine ver_hor_ad !$OMP PRIVATE(k) !$OMP DO do k=1,grd%km - grd%chl(:,:,k,l) = grd%chl(:,:,k,l) * grd%scy(:,:) + grd%chl(:,:,k,l) = grd%chl(:,:,k,l) * grd%scy(:,:,k) !laura enddo !$OMP END DO !$OMP END PARALLEL @@ -135,7 +135,7 @@ subroutine ver_hor_ad !$OMP PRIVATE(k) !$OMP DO do k=1,grd%km - grd%chl_ad(:,:,k,l) = grd%chl_ad(:,:,k,l) * grd%scy(:,:) + grd%chl_ad(:,:,k,l) = grd%chl_ad(:,:,k,l) * grd%scy(:,:,k) !laura enddo !$OMP END DO !$OMP END PARALLEL @@ -153,7 +153,7 @@ subroutine ver_hor_ad !$OMP PRIVATE(k) !$OMP DO do k=1,grd%km - grd%chl_ad(:,:,k,l) = grd%chl_ad(:,:,k,l) * grd%scx(:,:) + grd%chl_ad(:,:,k,l) = grd%chl_ad(:,:,k,l) * grd%scx(:,:,k) !laura enddo !$OMP END DO !$OMP END PARALLEL diff --git a/wrt_dia.f90 b/wrt_dia.f90 index 7fbeb5f..b48340a 100644 --- a/wrt_dia.f90 +++ b/wrt_dia.f90 @@ -62,7 +62,7 @@ subroutine wrt_dia do i=1,grd%im if (drv%argo .eq. 1) then if (grd%msk(i,j,k) .eq. 0) then - Dump_chl(i,j,k) = -1. + Dump_chl(i,j,k) = -1.!d+20 else Dump_chl(i,j,k) = REAL(grd%chl(i,j,k,l), 4) endif