From 91bc2cc18603796f4b2bb07873a06d14491fec40 Mon Sep 17 00:00:00 2001 From: Christopher Date: Thu, 4 Jul 2024 11:13:42 +0200 Subject: [PATCH 1/3] Initialize insig variables and add debug output --- src/tides.f90 | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/tides.f90 b/src/tides.f90 index d4bf9e9..e7e5c7d 100644 --- a/src/tides.f90 +++ b/src/tides.f90 @@ -479,6 +479,10 @@ subroutine getTidesAcceleration( & tanphi = tan(phi_gc) + insig1 = 0.d0 + insig2 = 0.d0 + insig3 = 0.d0 + do l = 2,lmax !** recursive computation of sin(ml) and cos(ml) @@ -512,6 +516,8 @@ subroutine getTidesAcceleration( & temp_t1 = -muEarth*oorabs3 temp_t2 = muEarth*oorabs + write(*,*) "temp_t = ", temp_t1, temp_t2, muEarth, oorabs3, oorabs + ! compute the radial partial derivative of the potential dudr = temp_t1 * insig1 @@ -521,6 +527,10 @@ subroutine getTidesAcceleration( & ! compute the longitudal partial derivative of the potential dudlambda = temp_t2 * insig3 + write(*,*) "du = ", dudr, dudphi, dudlambda + + write(*,*) "r_gcrf = ", r_gcrf + ! pre-compute terms which are used for the acceleration components r1r2 = r_gcrf(1)*r_gcrf(1) + r_gcrf(2)*r_gcrf(2) temp_t3 = 1.d0/r1r2 @@ -531,6 +541,8 @@ subroutine getTidesAcceleration( & temp2 = dudlambda*temp_t3 temp = dudr - temp_t4*dudphi + write(*,*) "temp = ", temp, temp2, temp_t5, dudr, dudphi + !========================================================================== ! ! Finally, compute the non-spherical perturbative accelerations in the GCRF @@ -546,8 +558,8 @@ subroutine getTidesAcceleration( & ! k-direction [km/s^2] accel(3) = dudr * r_gcrf(3) + temp_t5 * dudphi - !write(*,*) "accel = ", accel - !read(*,*) + write(*,*) "accel = ", accel, temp, temp2, temp_t5, dudr + read(*,*) !** done! if(isControlled()) then From fad6a323c5399aea4a0e4b7f00465405f76b05d5 Mon Sep 17 00:00:00 2001 From: Christopher Date: Thu, 4 Jul 2024 11:32:41 +0200 Subject: [PATCH 2/3] Adding more debug messages --- src/tides.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/tides.f90 b/src/tides.f90 index e7e5c7d..8add714 100644 --- a/src/tides.f90 +++ b/src/tides.f90 @@ -428,6 +428,8 @@ subroutine getTidesAcceleration( & end if + write (*,*) "dC, dS = ", dC, dS + !=============================================================================================== ! ! Compute required quantities required for the accelerations (or take from geoopotential...) @@ -440,6 +442,8 @@ subroutine getTidesAcceleration( & !** get radius, geocentric longitude and latitude call getRadiusLatLon(r_itrf, v_itrf, rabs, phi_gc, lambda) + write (*,*) "rabs, phi_gc, lambda = ", rabs, phi_gc, lambda + !** orbital radius rabs2 = rabs*rabs rabs3 = rabs2*rabs @@ -469,6 +473,8 @@ subroutine getTidesAcceleration( & end do + write (*,*) "lpsat = ", lpsat + ! determine partial derivatives of the disturbing potential costerm(0) = 1.d0 @@ -479,10 +485,6 @@ subroutine getTidesAcceleration( & tanphi = tan(phi_gc) - insig1 = 0.d0 - insig2 = 0.d0 - insig3 = 0.d0 - do l = 2,lmax !** recursive computation of sin(ml) and cos(ml) From 09272bfc7cb43b861e29e8e9703a953a7cf2f8f3 Mon Sep 17 00:00:00 2001 From: Christopher Date: Thu, 4 Jul 2024 12:06:12 +0200 Subject: [PATCH 3/3] Initialize dC and dS arrays and remove debug messages --- src/tides.f90 | 37 +++++++++++++------------------------ 1 file changed, 13 insertions(+), 24 deletions(-) diff --git a/src/tides.f90 b/src/tides.f90 index 8add714..df6c4c2 100644 --- a/src/tides.f90 +++ b/src/tides.f90 @@ -428,8 +428,6 @@ subroutine getTidesAcceleration( & end if - write (*,*) "dC, dS = ", dC, dS - !=============================================================================================== ! ! Compute required quantities required for the accelerations (or take from geoopotential...) @@ -442,8 +440,6 @@ subroutine getTidesAcceleration( & !** get radius, geocentric longitude and latitude call getRadiusLatLon(r_itrf, v_itrf, rabs, phi_gc, lambda) - write (*,*) "rabs, phi_gc, lambda = ", rabs, phi_gc, lambda - !** orbital radius rabs2 = rabs*rabs rabs3 = rabs2*rabs @@ -473,8 +469,6 @@ subroutine getTidesAcceleration( & end do - write (*,*) "lpsat = ", lpsat - ! determine partial derivatives of the disturbing potential costerm(0) = 1.d0 @@ -518,8 +512,6 @@ subroutine getTidesAcceleration( & temp_t1 = -muEarth*oorabs3 temp_t2 = muEarth*oorabs - write(*,*) "temp_t = ", temp_t1, temp_t2, muEarth, oorabs3, oorabs - ! compute the radial partial derivative of the potential dudr = temp_t1 * insig1 @@ -529,10 +521,6 @@ subroutine getTidesAcceleration( & ! compute the longitudal partial derivative of the potential dudlambda = temp_t2 * insig3 - write(*,*) "du = ", dudr, dudphi, dudlambda - - write(*,*) "r_gcrf = ", r_gcrf - ! pre-compute terms which are used for the acceleration components r1r2 = r_gcrf(1)*r_gcrf(1) + r_gcrf(2)*r_gcrf(2) temp_t3 = 1.d0/r1r2 @@ -543,8 +531,6 @@ subroutine getTidesAcceleration( & temp2 = dudlambda*temp_t3 temp = dudr - temp_t4*dudphi - write(*,*) "temp = ", temp, temp2, temp_t5, dudr, dudphi - !========================================================================== ! ! Finally, compute the non-spherical perturbative accelerations in the GCRF @@ -560,9 +546,6 @@ subroutine getTidesAcceleration( & ! k-direction [km/s^2] accel(3) = dudr * r_gcrf(3) + temp_t5 * dudphi - write(*,*) "accel = ", accel, temp, temp2, temp_t5, dudr - read(*,*) - !** done! if(isControlled()) then call checkOut(csubid) @@ -705,6 +688,11 @@ subroutine init_FES2004(this) real(dp) :: temp_dSm ! temporary retrograde S coefficient real(dp) :: temp_Doodson ! temporary Doodson number of the ocean tide + this%dC_p = 0.d0 + this%dS_p = 0.d0 + this%dC_m = 0.d0 + this%dS_m = 0.d0 + !read data from external file ich = openFile(trim(adjustl(this%dataPath))//cdelimit//trim(adjustl(this%fesDataFile)), SEQUENTIAL, IN_FORMATTED) do ind = 1, 59462 @@ -826,6 +814,7 @@ subroutine init_FES2004(this) end if end if end do + ich = closeFile(ich) !assign the right order of magnitude to harmonic coefficients this%dC_p = this%dC_p * 1.d-11 @@ -870,11 +859,11 @@ subroutine get_FES2004_corrections(this, time_mjd, lmax, reduction, dC, dS) real(dp), dimension(2:6,0:6), intent(inout) :: dC ! matrix with corrections on the C harmonic coefficients real(dp), dimension(2:6,0:6), intent(inout) :: dS ! matrix with corrections on the S harmonic coefficients - !get Delaunay arguments in radians for the current time + ! get Delaunay arguments in radians for the current time call get_Delaunay_arg(this, time_mjd, F_vect) F_vect = F_vect*deg2rad - !get Greenwich Mean Sidereal Time + ! get Greenwich Mean Sidereal Time theta_g = getGMST(time_mjd) do l = 2, lmax @@ -884,19 +873,19 @@ subroutine get_FES2004_corrections(this, time_mjd, lmax, reduction, dC, dS) else dm = 2 end if - !compute the corresponding unnormalization factor + ! compute the corresponding unnormalization factor fac = sqrt(factorial(l - m)*dm*(2.d0*l + 1.d0)/factorial(l + m)) do i = 1, 18 - !compute the argument for each tide constituent + !c ompute the argument for each tide constituent call get_Doodson_arg(this, this%doodson(i), ctheta_g, cl, cl_prime, cF, cD, cOmega) + theta_f = ctheta_g * theta_g + cl * F_vect(1) + cl_prime * F_vect(2) + cF * F_vect(3) + & cD * F_vect(4) + cOmega * F_vect(5) - - !get the produced gravity field corrections + ! get the produced gravity field corrections dC(l, m) = dC(l, m) + fac*((this%dC_p(i, l, m) + this%dC_m(i, l, m))*cos(theta_f) + & (this%dS_p(i, l, m) + this%dS_m(i, l, m))*sin(theta_f)) if (m == 0) then - dS(l, m) = 0 + dS(l, m) = 0.d0 else dS(l, m) = dS(l, m) + fac*((this%dS_p(i, l, m) - this%dS_m(i, l, m))*cos(theta_f) - & (this%dC_p(i, l, m) - this%dC_m(i, l, m))*sin(theta_f))