| 1 | subroutine s1_bar_func(xprime,yprime,zprime,s1_bar, |
|---|
| 2 | * s1_bar_x,s1_bar_y,s1_bar_z,g2s1b, |
|---|
| 3 | * Lz,s1_scale) |
|---|
| 4 | c |
|---|
| 5 | c *****input arguments xprime,yprime,zprime are DIMENSIONLESS as are all returned values |
|---|
| 6 | c (xprime,yprime,zprime): d'less spatial position |
|---|
| 7 | c Lz: vertical size of computational domain [m] |
|---|
| 8 | c s1_scale: characteristic range of s1 values, used to scale s1 |
|---|
| 9 | c units are dimensional, e.g. deg C for s1-->Temp. |
|---|
| 10 | |
|---|
| 11 | c *****output variables, these need to be specified in this routine |
|---|
| 12 | c s1_bar, s1_bar_x, s1_bar_y, s1_bar_z, g2s1b |
|---|
| 13 | c scalar concentration, and spatial gradients, g2s1b=grad2(s1_bar) [dimensional units] |
|---|
| 14 | |
|---|
| 15 | implicit none |
|---|
| 16 | real xprime,yprime,zprime |
|---|
| 17 | real Lz,s1_scale |
|---|
| 18 | real s1_bar,s1_bar_x,s1_bar_y,s1_bar_z,g2s1b |
|---|
| 19 | |
|---|
| 20 | c local variables: you may need some intermediate variables, |
|---|
| 21 | c these need to be explicitly declared here |
|---|
| 22 | real x,y,z,ans,ansp,ansm,delta_z,delta_x |
|---|
| 23 | |
|---|
| 24 | c it may be convenient to know the dimensional spatial location (x,y,z) in [m], so we compute it |
|---|
| 25 | x = xprime*Lz ! use dimensional value to evaluate formulae |
|---|
| 26 | y = yprime*Lz ! use dimensional value to evaluate formulae |
|---|
| 27 | z = zprime*Lz ! use dimensional value to evaluate formulae |
|---|
| 28 | |
|---|
| 29 | |
|---|
| 30 | c**** USER CUSTOMIZED CODE SEGMENT STARTS HERE ************************ |
|---|
| 31 | c scalars='TS' --> s1 = T' and s2 = S' |
|---|
| 32 | |
|---|
| 33 | delta_z = Lz/10000. ! [m] (offset for estimating derivs) |
|---|
| 34 | delta_x = 100. ! "" |
|---|
| 35 | |
|---|
| 36 | ! stay underwater |
|---|
| 37 | if(z >= Lz-delta_z) z=Lz-delta_z |
|---|
| 38 | |
|---|
| 39 | call interp2d(x,z,'T',ans) |
|---|
| 40 | call interp2d(x,z+delta_z,'T',ansp) |
|---|
| 41 | call interp2d(x,z-delta_z,'T',ansm) |
|---|
| 42 | cc print*,'after call interp2d on z' |
|---|
| 43 | |
|---|
| 44 | s1_bar = ans ! [deg C] |
|---|
| 45 | s1_bar_z = (ansp-ansm)/(2.*delta_z) ! [deg C/m] |
|---|
| 46 | g2s1b = (ansp -2.*ans + ansm)/(delta_z**2) ! [deg C/m2] |
|---|
| 47 | |
|---|
| 48 | |
|---|
| 49 | call interp2d(x+delta_x,z,'T',ansp) |
|---|
| 50 | call interp2d(x-delta_x,z,'T',ansm) |
|---|
| 51 | cc print*,'after call interp2d on x' |
|---|
| 52 | |
|---|
| 53 | |
|---|
| 54 | s1_bar_x = (ansp-ansm)/(2.*delta_x) ! [deg C/m] |
|---|
| 55 | |
|---|
| 56 | |
|---|
| 57 | s1_bar_y = 0.0 ! [deg C/m] |
|---|
| 58 | |
|---|
| 59 | |
|---|
| 60 | c**** USER DEFINED CODE SEGMENT ENDS HERE ******************************** |
|---|
| 61 | |
|---|
| 62 | |
|---|
| 63 | |
|---|
| 64 | |
|---|
| 65 | c*** LEAVE THIS SECTION AS IS, SCALING MUST BE CONSISTENT******************** |
|---|
| 66 | c*** WITH THE REST OF THE ALGORITHM **************************************** |
|---|
| 67 | s1_bar=s1_bar/(s1_scale) |
|---|
| 68 | s1_bar_x = s1_bar_x/(s1_scale/Lz) |
|---|
| 69 | s1_bar_y = s1_bar_y/(s1_scale/Lz) |
|---|
| 70 | s1_bar_z = s1_bar_z/(s1_scale/Lz) |
|---|
| 71 | g2s1b = g2s1b/(s1_scale/Lz**2) |
|---|
| 72 | c***************************************************************************** |
|---|
| 73 | c***************************************************************************** |
|---|
| 74 | |
|---|
| 75 | c print*,'s1_bar interpoled ...' |
|---|
| 76 | return |
|---|
| 77 | end |
|---|
| 78 | |
|---|
| 79 | c It may seem odd to have the inputs in dimensionless form |
|---|
| 80 | c when it is usually convenient for the user to define the |
|---|
| 81 | c ambient quantities in dimensional form. Also, we have a |
|---|
| 82 | c **do not touch** section of code that converts the user |
|---|
| 83 | c specified dimensional quantities back to dimensionless form. |
|---|
| 84 | c I chose to do it this way because this function gets called |
|---|
| 85 | c from so many places in the code, it would be a pain and |
|---|
| 86 | c a potential source of typos to do the scalings just before |
|---|
| 87 | c and just after each call. This is not the case with |
|---|
| 88 | c user_defined_ics and user_defined_forcing. These routines |
|---|
| 89 | c get called in one location only and so it makes sense to |
|---|
| 90 | c make the scaling "invisible" to the user. |
|---|
