anchor.diff:  
       ! magnitude term of gaussian at x0
       deltav0m = wind_speed * (                                         &
  -                (1.0_dp - sqrt(1.0_dp - Ct *                     &
  -                cos(yaw) / (8.0_dp * sigmay0 * sigmaz0 /           &
  +                (1.0_dp - sqrt(1.0_dp - Ct *                          &
  +                cos(yaw) / (8.0_dp * sigmay0 * sigmaz0 /              &
                                               (rotor_diameter ** 2)))))
                                               
       ! linearized gaussian magnitude term for near wake
  -    deltav = (((deltav0m - deltavs)/x0) * x + deltavs) *  &
  -        exp(-0.5_dp * (deltay / sigmay) ** 2) *                   &
  +    deltav = (((deltav0m - deltavs)/x0) * x + deltavs) *              &
  +        exp(-0.5_dp * (deltay / sigmay) ** 2) *                       &
           exp(-0.5_dp * (deltaz / sigmaz) ** 2)
       
   end subroutine deltav_near_wake_lin_func
  -
  -
  -!        Generated by TAPENADE     (INRIA, Ecuador team)
  -!  Tapenade 3.12 (r6213) - 13 Oct 2016 10:54
  -!
  -!  Differentiation of porteagel_analyze in reverse (adjoint) mode:
  -!   gradient     of useful results: wtvelocity
  -!   with respect to varying inputs: rotordiameter turbinez turbinexw
  -!                wtvelocity turbineyw yawdeg ct
  -!   RW status of diff variables: rotordiameter:out turbinez:out
  -!                turbinexw:out wtvelocity:in-zero turbineyw:out
  -!                yawdeg:out ct:out
  -SUBROUTINE PORTEAGEL_ANALYZE_BV(nturbines, turbinexw, turbinexwb, &
  -& turbineyw, turbineywb, turbinez, turbinezb, rotordiameter, &
  -& rotordiameterb, ct, ctb, wind_speed, yawdeg, yawdegb, ky, kz, alpha, &
  -& beta, i, wtvelocity, wtvelocityb, nbdirs)
  -!  Hint: nbdirs should be the maximum number of differentiation directions
  -  IMPLICIT NONE
  -!print *, "fortran"
  -! define precision to be the standard for a double precision ! on local system
  -  INTEGER, PARAMETER :: dp=KIND(0.d0)
  -! in
  -  INTEGER, INTENT(IN) :: nturbines
  -  REAL(dp), DIMENSION(nturbines), INTENT(IN) :: turbinexw, turbineyw, &
  -& turbinez
  -  REAL(dp), DIMENSION(nbdirs, nturbines), intent(out) :: turbinexwb, turbineywb, &
  -& turbinezb
  -  REAL(dp), DIMENSION(nturbines), INTENT(IN) :: rotordiameter, yawdeg
  -  REAL(dp), DIMENSION(nbdirs, nturbines), intent(out) :: rotordiameterb, yawdegb
  -  REAL(dp), DIMENSION(nturbines), INTENT(IN) :: ct
  -  REAL(dp), DIMENSION(nbdirs, nturbines), intent(out) :: ctb
  -  REAL(dp), INTENT(IN) :: ky, kz, alpha, beta, i, wind_speed
  -! local (General)
  -  REAL(dp), DIMENSION(nturbines) :: yaw
  -  REAL(dp), DIMENSION(nbdirs, nturbines) :: yawb
  -  REAL(dp) :: x0, deltax0, deltay, theta_c_0, sigmay, sigmaz, &
  -& wake_offset
  -  REAL(dp), DIMENSION(nbdirs) :: x0b, deltax0b, deltayb, theta_c_0b, &
  -& sigmayb, sigmazb, wake_offsetb
  -  REAL(dp) :: deltav
  -  REAL(dp), DIMENSION(nbdirs) :: deltavb
  -  INTEGER :: turb, turbi
  -  REAL(dp), PARAMETER :: pi=3.141592653589793_dp
  -! model out
  -  REAL(dp), DIMENSION(nturbines), intent(out) :: wtvelocity
  -  REAL(dp), DIMENSION(nbdirs, nturbines) :: wtvelocityb
  -  INTRINSIC COS, ATAN, MAX, SQRT, LOG
  -  INTRINSIC KIND
  -  INTRINSIC EXP
  -  INTEGER :: nd
  -  REAL(dp) :: temp
  -  REAL(dp) :: temp0
  -  REAL(dp) :: temp1
  -  REAL(dp) :: temp2
  -  REAL(dp) :: temp3
  -  REAL(dp) :: temp4
  -  REAL(dp) :: temp5
  -  REAL(dp) :: temp6
  -  REAL(dp) :: temp7
  -  REAL(dp) :: temp8
  -  REAL(dp), DIMENSION(nbdirs) :: tempb
  -  REAL(dp), DIMENSION(nbdirs) :: tempb0
  -  REAL(dp), DIMENSION(nbdirs) :: tempb1
  -  REAL(dp) :: temp9
  -  REAL(dp) :: temp10
  -  REAL(dp) :: temp11
  -  REAL(dp) :: temp12
  -  REAL(dp) :: temp13
  -  REAL(dp) :: temp14
  -  REAL(dp) :: temp15
  -  REAL(dp) :: temp16
  -  REAL(dp) :: temp17
  -  REAL(dp) :: temp18
  -  REAL(dp) :: temp19
  -  REAL(dp) :: temp20
  -  REAL(dp) :: temp21
  -  REAL(dp) :: temp22
  -  REAL(dp) :: temp23
  -  REAL(dp) :: temp24
  -  REAL(dp) :: temp25
  -  REAL(dp) :: temp26
  -  REAL(dp) :: temp27
  -  REAL(dp) :: temp28
  -  REAL(dp) :: temp29
  -  REAL(dp) :: temp30
  -  REAL(dp) :: temp31
  -  REAL(dp) :: temp32
  -  REAL(dp) :: temp33
  -  REAL(dp) :: temp34
  -  REAL(dp) :: temp35
  -  REAL(dp) :: temp36
  -  REAL(dp) :: temp37
  -  REAL(dp) :: temp38
  -  REAL(dp) :: temp39
  -  REAL(dp) :: temp40
  -  REAL(dp) :: temp41
  -  REAL(dp) :: temp42
  -  REAL(dp) :: temp43
  -  REAL(dp) :: temp44
  -  REAL(dp) :: temp45
  -  REAL(dp) :: temp46
  -  REAL(dp) :: temp47
  -  REAL(dp), DIMENSION(nbdirs) :: tempb2
  -  REAL(dp), DIMENSION(nbdirs) :: tempb3
  -  REAL(dp), DIMENSION(nbdirs) :: tempb4
  -  REAL(dp), DIMENSION(nbdirs) :: tempb5
  -  REAL(dp), DIMENSION(nbdirs) :: tempb6
  -  REAL(dp), DIMENSION(nbdirs) :: tempb7
  -  REAL(dp), DIMENSION(nbdirs) :: tempb8
  -  REAL(dp), DIMENSION(nbdirs) :: tempb9
  -  REAL(dp), DIMENSION(nbdirs) :: tempb10
  -  REAL(dp), DIMENSION(nbdirs) :: tempb11
  -  REAL(dp), DIMENSION(nbdirs) :: tempb12
  -  REAL(dp), DIMENSION(nbdirs) :: tempb13
  -  REAL(dp), DIMENSION(nbdirs) :: tempb14
  -  REAL(dp), DIMENSION(nbdirs) :: tempb15
  -  REAL(dp), DIMENSION(nbdirs) :: tempb16
  -  REAL(dp), DIMENSION(nbdirs) :: tempb17
  -  REAL(dp), DIMENSION(nbdirs) :: tempb18
  -  REAL(dp), DIMENSION(nbdirs) :: tempb19
  -  REAL(dp), DIMENSION(nbdirs) :: tempb20
  -  REAL(dp), DIMENSION(nbdirs) :: tempb21
  -  INTEGER :: branch
  -  INTEGER :: nbdirs
  -  yaw = -(yawdeg*pi/180.0_dp)
  -  DO turb=1,nturbines
  -    CALL PUSHREAL4ARRAY(x0, dp/4)
  -    x0 = rotordiameter(turb)*(COS(yaw(turb))*(1.0_dp+SQRT(1.0_dp-ct(turb&
  -&     )))/(SQRT(2.0_dp)*(alpha*i+beta*(1.0_dp-SQRT(1.0_dp-ct(turb))))))
  -    CALL PUSHREAL4ARRAY(theta_c_0, dp/4)
  -    theta_c_0 = 0.3_dp*yaw(turb)*(1.0_dp-SQRT(1.0_dp-ct(turb)*COS(yaw(&
  -&     turb))))/COS(yaw(turb))
  -! at turbineX-locations
  -    DO turbi=1,nturbines
  -      CALL PUSHREAL4ARRAY(deltax0, dp/4)
  -      deltax0 = turbinexw(turbi) - (turbinexw(turb)+x0)
  -      IF (deltax0 .GT. 0.0_dp) THEN
  -        CALL PUSHREAL4ARRAY(sigmay, dp/4)
  -        sigmay = rotordiameter(turb)*(ky*deltax0/rotordiameter(turb)+COS&
  -&         (yaw(turb))/SQRT(8.0_dp))
  -        CALL PUSHREAL4ARRAY(sigmaz, dp/4)
  -        sigmaz = rotordiameter(turb)*(kz*deltax0/rotordiameter(turb)+&
  -&         1.0_dp/SQRT(8.0_dp))
  -        wake_offset = rotordiameter(turb)*(theta_c_0*x0/rotordiameter(&
  -&         turb)+theta_c_0/14.7_dp*SQRT(COS(yaw(turb))/(ky*kz*ct(turb)))*&
  -&         (2.9_dp+1.3_dp*SQRT(1.0_dp-ct(turb))-ct(turb))*LOG((1.6_dp+&
  -&         SQRT(ct(turb)))*(1.6_dp*SQRT(8.0_dp*sigmay*sigmaz/(COS(yaw(&
  -&         turb))*rotordiameter(turb)**2))-SQRT(ct(turb)))/((1.6_dp-SQRT(&
  -&         ct(turb)))*(1.6_dp*SQRT(8.0_dp*sigmay*sigmaz/(COS(yaw(turb))*&
  -&         rotordiameter(turb)**2))+SQRT(ct(turb))))))
  -        CALL PUSHREAL4ARRAY(deltay, dp/4)
  -        deltay = turbineyw(turbi) - (turbineyw(turb)+wake_offset)
  -        CALL PUSHCONTROL1B(1)
  -      ELSE
  -        CALL PUSHCONTROL1B(0)
  -      END IF
  -    END DO
  -  END DO
  -  DO nd=1,nbdirs
  -    rotordiameterb(nd, :) = 0.0
  -    turbinezb(nd, :) = 0.0
  -    turbinexwb(nd, :) = 0.0
  -    turbineywb(nd, :) = 0.0
  -    ctb(nd, :) = 0.0
  -    yawb(nd, :) = 0.0
  -  END DO
  -  DO turb=nturbines,1,-1
  -    DO nd=1,nbdirs
  -      theta_c_0b(nd) = 0.0
  -      x0b(nd) = 0.0
  -    END DO
  -    DO turbi=nturbines,1,-1
  -      CALL POPCONTROL1B(branch)
  -      IF (branch .EQ. 0) THEN
  -        DO nd=1,nbdirs
  -          deltax0b(nd) = 0.0
  -        END DO
  -      ELSE
  -        temp47 = sigmaz**2
  -        temp40 = turbinez(turbi) - turbinez(turb)
  -        temp39 = temp40**2/temp47
  -        temp46 = EXP(-(0.5_dp*temp39))
  -        temp45 = sigmay**2
  -        temp38 = deltay**2/temp45
  -        temp44 = EXP(-(0.5_dp*temp38))
  -        temp37 = 8.0_dp*sigmay*sigmaz
  -        temp43 = ct(turb)*rotordiameter(turb)**2
  -        temp36 = temp43/temp37
  -        temp42 = COS(yaw(turb))
  -        temp41 = SQRT(-(temp42*temp36) + 1.0_dp)
  -        temp35 = SQRT(ct(turb))
  -        temp34 = COS(yaw(turb))
  -        temp20 = temp34*rotordiameter(turb)**2
  -        temp19 = sigmay*sigmaz/temp20
  -        temp33 = SQRT(8.0_dp*temp19)
  -        temp32 = SQRT(ct(turb))
  -        temp18 = (-temp32+1.6_dp)*(1.6_dp*temp33+temp35)
  -        temp31 = SQRT(ct(turb))
  -        temp30 = COS(yaw(turb))
  -        temp17 = temp30*rotordiameter(turb)**2
  -        temp16 = sigmay*sigmaz/temp17
  -        temp29 = SQRT(8.0_dp*temp16)
  -        temp28 = 1.6_dp*temp29 - temp31
  -        temp27 = SQRT(ct(turb))
  -        temp15 = (temp27+1.6_dp)*temp28/temp18
  -        temp26 = LOG(temp15)
  -        temp25 = SQRT(-ct(turb) + 1.0_dp)
  -        temp14 = 1.3_dp*temp25 - ct(turb) + 2.9_dp
  -        temp24 = theta_c_0*temp14/14.7_dp
  -        temp23 = ky*kz*ct(turb)
  -        temp22 = COS(yaw(turb))
  -        temp13 = temp22/temp23
  -        temp21 = SQRT(temp13)
  -        temp12 = theta_c_0*x0/rotordiameter(turb)
  -        deltax0 = turbinexw(turbi) - (turbinexw(turb)+x0)
  -        temp11 = deltax0/rotordiameter(turb)
  -        temp10 = SQRT(8.0_dp)
  -        temp9 = deltax0/rotordiameter(turb)
  -        DO nd=1,nbdirs
  -          deltavb(nd) = -wtvelocityb(nd, turbi)
  -          IF (1.0_dp - temp42*temp36 .EQ. 0.0) THEN
  -            tempb2(nd) = 0.0
  -          ELSE
  -            tempb2(nd) = -(temp44*temp46*wind_speed*deltavb(nd)/(2.0*&
  -&             temp41))
  -          END IF
  -          tempb3(nd) = -(temp42*tempb2(nd)/temp37)
  -          tempb4(nd) = -(temp36*tempb3(nd))
  -          tempb5(nd) = wind_speed*(1.0_dp-temp41)*deltavb(nd)
  -          tempb6(nd) = -(0.5_dp*EXP(-(0.5_dp*temp38))*temp46*tempb5(nd)/&
  -&           temp45)
  -          tempb7(nd) = -(0.5_dp*EXP(-(0.5_dp*temp39))*temp44*tempb5(nd)/&
  -&           temp47)
  -          tempb8(nd) = 2*temp40*tempb7(nd)
  -          ctb(nd, turb) = ctb(nd, turb) + rotordiameter(turb)**2*tempb3(&
  -&           nd)
  -          deltayb(nd) = 2*deltay*tempb6(nd)
  -          turbinezb(nd, turbi) = turbinezb(nd, turbi) + tempb8(nd)
  -          turbinezb(nd, turb) = turbinezb(nd, turb) - tempb8(nd)
  -          turbineywb(nd, turbi) = turbineywb(nd, turbi) + deltayb(nd)
  -          turbineywb(nd, turb) = turbineywb(nd, turb) - deltayb(nd)
  -          wake_offsetb(nd) = -deltayb(nd)
  -          tempb17(nd) = rotordiameter(turb)*wake_offsetb(nd)
  -          tempb12(nd) = tempb17(nd)/rotordiameter(turb)
  -          tempb18(nd) = temp26*tempb17(nd)
  -          IF (temp13 .EQ. 0.0) THEN
  -            tempb9(nd) = 0.0
  -          ELSE
  -            tempb9(nd) = temp24*tempb18(nd)/(2.0*temp21*temp23)
  -          END IF
  -          tempb19(nd) = temp21*theta_c_0*tempb18(nd)/14.7_dp
  -          tempb20(nd) = temp21*temp24*tempb17(nd)/(temp15*temp18)
  -          IF (8.0_dp*temp16 .EQ. 0.0) THEN
  -            tempb15(nd) = 0.0
  -          ELSE
  -            tempb15(nd) = 8.0_dp*1.6_dp*(temp27+1.6_dp)*tempb20(nd)/(2.0&
  -&             *temp29*temp17)
  -          END IF
  -          tempb10(nd) = -(temp16*tempb15(nd))
  -          tempb21(nd) = -(temp15*tempb20(nd))
  -          IF (8.0_dp*temp19 .EQ. 0.0) THEN
  -            tempb16(nd) = 0.0
  -          ELSE
  -            tempb16(nd) = 8.0_dp*1.6_dp*(1.6_dp-temp32)*tempb21(nd)/(2.0&
  -&             *temp33*temp20)
  -          END IF
  -          sigmayb(nd) = sigmaz*tempb15(nd) + sigmaz*tempb16(nd) - temp38&
  -&           *2*sigmay*tempb6(nd) + sigmaz*8.0_dp*tempb4(nd)
  -          sigmazb(nd) = sigmay*tempb15(nd) + sigmay*tempb16(nd) - temp39&
  -&           *2*sigmaz*tempb7(nd) + 8.0_dp*sigmay*tempb4(nd)
  -          tempb11(nd) = -(temp19*tempb16(nd))
  -          yawb(nd, turb) = yawb(nd, turb) + temp36*SIN(yaw(turb))*tempb2&
  -&           (nd) - rotordiameter(turb)**2*SIN(yaw(turb))*tempb10(nd) - &
  -&           rotordiameter(turb)**2*SIN(yaw(turb))*tempb11(nd) - &
  -&           rotordiameter(turb)*SIN(yaw(turb))*sigmayb(nd)/temp10 - SIN(&
  -&           yaw(turb))*tempb9(nd)
  -          theta_c_0b(nd) = theta_c_0b(nd) + temp14*temp21*tempb18(nd)/&
  -&           14.7_dp + x0*tempb12(nd)
  -          x0b(nd) = x0b(nd) + theta_c_0*tempb12(nd)
  -          IF (1.0_dp - ct(turb) .EQ. 0.0) THEN
  -            ctb(nd, turb) = ctb(nd, turb) - tempb19(nd) - temp13*ky*kz*&
  -&             tempb9(nd)
  -          ELSE
  -            ctb(nd, turb) = ctb(nd, turb) + ((1.6_dp-temp32)/(2.0*temp35&
  -&             )-(1.6_dp*temp33+temp35)/(2.0*temp32))*tempb21(nd) + (&
  -&             temp28/(2.0*temp27)-(temp27+1.6_dp)/(2.0*temp31))*tempb20(&
  -&             nd) + ((-1.0)-1.3_dp/(2.0*temp25))*tempb19(nd) - temp13*ky&
  -&             *kz*tempb9(nd)
  -          END IF
  -          tempb14(nd) = kz*sigmazb(nd)
  -          tempb13(nd) = ky*sigmayb(nd)
  -          rotordiameterb(nd, turb) = rotordiameterb(nd, turb) + (temp12+&
  -&           temp21*temp24*temp26)*wake_offsetb(nd) - temp12*tempb12(nd) &
  -&           + temp30*2*rotordiameter(turb)*tempb10(nd) + temp34*2*&
  -&           rotordiameter(turb)*tempb11(nd) + (ky*temp9+COS(yaw(turb))/&
  -&           temp10)*sigmayb(nd) - temp9*tempb13(nd) - temp11*tempb14(nd)&
  -&           + (1.0/SQRT(8.0_dp)+kz*temp11)*sigmazb(nd) + ct(turb)*2*&
  -&           rotordiameter(turb)*tempb3(nd)
  -          deltax0b(nd) = tempb13(nd) + tempb14(nd)
  -        END DO
  -        CALL POPREAL4ARRAY(deltay, dp/4)
  -        CALL POPREAL4ARRAY(sigmaz, dp/4)
  -        CALL POPREAL4ARRAY(sigmay, dp/4)
  -      END IF
  -      CALL POPREAL4ARRAY(deltax0, dp/4)
  -      DO nd=1,nbdirs
  -        turbinexwb(nd, turbi) = turbinexwb(nd, turbi) + deltax0b(nd)
  -        turbinexwb(nd, turb) = turbinexwb(nd, turb) - deltax0b(nd)
  -        x0b(nd) = x0b(nd) - deltax0b(nd)
  -      END DO
  -    END DO
  -    CALL POPREAL4ARRAY(theta_c_0, dp/4)
  -    temp8 = COS(yaw(turb))
  -    temp5 = yaw(turb)/temp8
  -    temp7 = COS(yaw(turb))
  -    temp4 = -(ct(turb)*temp7) + 1.0_dp
  -    temp6 = SQRT(temp4)
  -    CALL POPREAL4ARRAY(x0, dp/4)
  -    temp3 = SQRT(-ct(turb) + 1.0_dp)
  -    temp2 = SQRT(2.0_dp)
  -    temp1 = temp2*(alpha*i+beta*(-temp3+1.0_dp))
  -    temp0 = SQRT(-ct(turb) + 1.0_dp)
  -    temp = COS(yaw(turb))
  -    DO nd=1,nbdirs
  -      tempb1(nd) = x0b(nd)/temp1
  -      IF (temp4 .EQ. 0.0) THEN
  -        tempb(nd) = 0.0
  -      ELSE
  -        tempb(nd) = -(temp5*0.3_dp*theta_c_0b(nd)/(2.0*temp6))
  -      END IF
  -      tempb0(nd) = (1.0_dp-temp6)*0.3_dp*theta_c_0b(nd)/temp8
  -      ctb(nd, turb) = ctb(nd, turb) - temp7*tempb(nd)
  -      yawb(nd, turb) = yawb(nd, turb) + (temp5*SIN(yaw(turb))+1.0)*&
  -&       tempb0(nd) - rotordiameter(turb)*(temp0+1.0_dp)*SIN(yaw(turb))*&
  -&       tempb1(nd) + ct(turb)*SIN(yaw(turb))*tempb(nd)
  -      rotordiameterb(nd, turb) = rotordiameterb(nd, turb) + (temp0+&
  -&       1.0_dp)*temp*tempb1(nd)
  -      IF (.NOT.1.0_dp - ct(turb) .EQ. 0.0) ctb(nd, turb) = ctb(nd, turb)&
  -&         + (-(beta*temp2*rotordiameter(turb)*temp*(temp0+1.0_dp)/(2.0*&
  -&         temp3*temp1))-rotordiameter(turb)*temp/(2.0*temp0))*tempb1(nd)
  -    END DO
  -  END DO
  -  DO nd=1,nbdirs
  -    yawdegb(nd, :) = 0.0
  -    yawdegb(nd, :) = -(pi*yawb(nd, :)/180.0_dp)
  -  END DO
  -  DO nd=1,nbdirs
  -    wtvelocityb(nd, :) = 0.0
  -  END DO
  -END SUBROUTINE PORTEAGEL_ANALYZE_BV
  
positive.code: 
  end subroutine porteagel_analyze
  
  
  subroutine porteagel_visualize(nTurbines, nSamples, turbineXw, turbineYw, turbineZ, &
                               rotorDiameter, Ct, wind_speed, &
                               yawDeg, ky, kz, alpha, beta, I, velX, velY, velZ, wsArray)
  
negative.code: ! implementation of the Bastankhah and Porte Agel (BPA) wake model for analysis
  subroutine porteagel_analyze(nTurbines, turbineXw, turbineYw, turbineZ, &
                               rotorDiameter, Ct, wind_speed, &
                               yawDeg, ky, kz, alpha, beta, I, wtVelocity)
  