Initialize dC and dS arrays and remove debug messages

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))