File indexing completed on 2023-09-03 05:10:16 UTC
view on githubraw file Latest commit 74487008 on 2023-09-03 01:50:18 UTC
6d54cf9ca1 Ed H*0001 #include "EXF_OPTIONS.h"
3a255f48df Gael*0002 #ifdef ALLOW_AUTODIFF
0003 # include "AUTODIFF_OPTIONS.h"
0004 #endif
bdec91d862 Patr*0005
d216b1690c Mart*0006
0007
0008
0009 SUBROUTINE EXF_BULKFORMULAE( myTime, myIter, myThid )
0010
0011
0012
0013
0014
0015
e1fb02e8f0 Jean*0016
d216b1690c Mart*0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141 IMPLICIT NONE
0142
bdec91d862 Patr*0143 #include "EEPARAMS.h"
0144 #include "SIZE.h"
0145 #include "PARAMS.h"
0146 #include "DYNVARS.h"
d216b1690c Mart*0147 #include "GRID.h"
bdec91d862 Patr*0148
082e18c36c Jean*0149 #include "EXF_PARAM.h"
0150 #include "EXF_FIELDS.h"
0151 #include "EXF_CONSTANTS.h"
bdec91d862 Patr*0152
0153 #ifdef ALLOW_AUTODIFF_TAMC
0154 #include "tamc.h"
0155 #endif
0156
d216b1690c Mart*0157
0158
0159
0160
0161
0162 _RL myTime
0163 INTEGER myIter
0164 INTEGER myThid
0165
0166
bdec91d862 Patr*0167
7109a141b2 Patr*0168 #ifdef ALLOW_BULKFORMULAE
423768d890 Jean*0169 #ifdef ALLOW_ATM_TEMP
d216b1690c Mart*0170
7109a141b2 Patr*0171
d216b1690c Mart*0172
0173
0174
0175
0176 INTEGER i,j,bi,bj
bdec91d862 Patr*0177
20cf30c902 Patr*0178 _RL czol
d216b1690c Mart*0179 _RL Tsf
0180 _RL wsm
0181 _RL t0
0182
0183 _RL tstar (1:sNx,1:sNy)
0184 _RL qstar (1:sNx,1:sNy)
0185 _RL ustar (1:sNx,1:sNy)
0186 _RL tau (1:sNx,1:sNy)
0187 _RL rdn (1:sNx,1:sNy)
0188 _RL rd (1:sNx,1:sNy)
0189 _RL delq (1:sNx,1:sNy)
0b902863ee Mart*0190 _RL deltap(1:sNx,1:sNy)
8053a926ae Jean*0191 #ifdef EXF_CALC_ATMRHO
0192 _RL atmrho_loc(1:sNx,1:sNy)
0193 #endif
0194
e1fb02e8f0 Jean*0195 #ifdef ALLOW_BULK_LARGEYEAGER04
d216b1690c Mart*0196 _RL dzTmp
e1fb02e8f0 Jean*0197 #endif
d216b1690c Mart*0198 _RL SSTtmp
0199 _RL ssq
0200 _RL re
0201 _RL rh
0202 _RL ren, rhn
0203 _RL usn, usm
0204 _RL stable
0205 _RL huol
0206 _RL htol
0207 _RL hqol
0208 _RL x
0209 _RL xsq
0210 _RL psimh
0211 _RL psixh
0212 _RL zwln
0213 _RL ztln
0214 _RL tmpbulk
0215 _RL recip_rhoConstFresh
0320e25227 Mart*0216 INTEGER ks, kl
d216b1690c Mart*0217 INTEGER iter
0218
0219 LOGICAL solve4Stress
e1fb02e8f0 Jean*0220 _RL windStress
d216b1690c Mart*0221
bdec91d862 Patr*0222 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0223 INTEGER ikey, ikey_1, ikey_2
bdec91d862 Patr*0224 #endif
0225
e1fb02e8f0 Jean*0226
0227
423768d890 Jean*0228 IF ( useAtmWind ) THEN
358649780a Gael*0229 solve4Stress = .TRUE.
423768d890 Jean*0230 ELSE
d216b1690c Mart*0231 #ifdef ALLOW_BULK_LARGEYEAGER04
358649780a Gael*0232 solve4Stress = wspeedfile .NE. ' '
d216b1690c Mart*0233 #else
358649780a Gael*0234 solve4Stress = .FALSE.
d216b1690c Mart*0235 #endif
423768d890 Jean*0236 ENDIF
bdec91d862 Patr*0237
0320e25227 Mart*0238 IF ( usingPCoords ) THEN
0239 ks = Nr
0240 kl = Nr-1
0241 ELSE
0242 ks = 1
0243 kl = 2
0244 ENDIF
0245
0246
666d41f610 Mart*0247 zwln = LOG(hu/zref)
0248 ztln = LOG(ht/zref)
20cf30c902 Patr*0249 czol = hu*karman*gravity_mks
d216b1690c Mart*0250
666d41f610 Mart*0251 recip_rhoConstFresh = 1. _d 0/rhoConstFresh
e1fb02e8f0 Jean*0252
8053a926ae Jean*0253
bdec91d862 Patr*0254 #ifdef ALLOW_AUTODIFF_TAMC
0255
0256
0257 #endif
d216b1690c Mart*0258 DO bj = myByLo(myThid),myByHi(myThid)
bdec91d862 Patr*0259 #ifdef ALLOW_AUTODIFF_TAMC
0260
0261
0262 #endif
d216b1690c Mart*0263 DO bi = myBxLo(myThid),myBxHi(myThid)
0320e25227 Mart*0264
20dee61641 Mart*0265 #ifdef ALLOW_AUTODIFF_TAMC
0266 ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
0267 #endif
d216b1690c Mart*0268 DO j = 1,sNy
0269 DO i = 1,sNx
bdec91d862 Patr*0270
0271 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0272 ikey_1 = i + (j-1)*sNx + (ikey-1)*sNx*sNy
bdec91d862 Patr*0273 #endif
0274
423768d890 Jean*0275 #ifdef ALLOW_AUTODIFF
f1644eaf0d Patr*0276 deltap(i,j) = 0. _d 0
0277 delq(i,j) = 0. _d 0
0278 #endif
d216b1690c Mart*0279
0280
0281 IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
0282
0320e25227 Mart*0283 Tsf = theta(i,j,ks,bi,bj) + cen2kel
c2c9eed210 Jean*0284 IF ( Nr.GE.2 .AND. sstExtrapol.GT.0. _d 0 ) THEN
d216b1690c Mart*0285 SSTtmp = sstExtrapol
0320e25227 Mart*0286 & *( theta(i,j,ks,bi,bj)-theta(i,j,kl,bi,bj) )
0287 & * maskC(i,j,kl,bi,bj)
d216b1690c Mart*0288 Tsf = Tsf + MAX( SSTtmp, 0. _d 0 )
0289 ENDIF
8053a926ae Jean*0290
d216b1690c Mart*0291 tmpbulk = cvapor_fac*exp(-cvapor_exp/Tsf)
8053a926ae Jean*0292 #ifdef EXF_CALC_ATMRHO
0293 atmrho_loc(i,j) = apressure(i,j,bi,bj) /
0294 & (287.04 _d 0*atemp(i,j,bi,bj)
0295 & *(1. _d 0 + humid_fac*aqh(i,j,bi,bj)))
0296 ssq = saltsat*tmpbulk/atmrho_loc(i,j)
0297 #else
d216b1690c Mart*0298 ssq = saltsat*tmpbulk/atmrho
8053a926ae Jean*0299 #endif
d216b1690c Mart*0300 deltap(i,j) = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
0b902863ee Mart*0301 delq(i,j) = aqh(i,j,bi,bj) - ssq
d216b1690c Mart*0302
0303 IF ( noNegativeEvap ) delq(i,j) = MIN( 0. _d 0, delq(i,j) )
0304
0305
0306
0307 stable = exf_half + SIGN(exf_half, deltap(i,j))
0308
0309 IF ( solve4Stress ) THEN
bdec91d862 Patr*0310 #ifdef ALLOW_AUTODIFF_TAMC
3752238fd8 Patr*0311
d216b1690c Mart*0312 #endif /* ALLOW_AUTODIFF_TAMC */
0313
0314 wsm = sh(i,j,bi,bj)
cc60455fbb Mart*0315 #ifdef ALLOW_DRAG_LARGEYEAGER09
0316
0317 tmpbulk = cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
0318 & + cdrag_8 * wsm**6
0319 tmpbulk = exf_scal_BulkCdn * (
0320 & ( halfRL - SIGN(halfRL, wsm-umax) )*tmpbulk
0321 & + ( halfRL + SIGN(halfRL, wsm-umax) )*cdragMax
0322 & )
0323 #else
d216b1690c Mart*0324 tmpbulk = exf_scal_BulkCdn *
0325 & ( cdrag_1/wsm + cdrag_2 + cdrag_3*wsm )
cc60455fbb Mart*0326 #endif
d216b1690c Mart*0327 rdn(i,j) = SQRT(tmpbulk)
0328 ustar(i,j) = rdn(i,j)*wsm
0329 ELSE
ddea17ebdb Gael*0330 rdn(i,j) = 0. _d 0
d216b1690c Mart*0331 windStress = wStress(i,j,bi,bj)
8053a926ae Jean*0332 #ifdef EXF_CALC_ATMRHO
423768d890 Jean*0333 ustar(i,j) = SQRT(windStress/atmrho_loc(i,j))
0334 tau(i,j) = SQRT(windStress*atmrho_loc(i,j))
8053a926ae Jean*0335 #else
423768d890 Jean*0336 ustar(i,j) = SQRT(windStress/atmrho)
d216b1690c Mart*0337
423768d890 Jean*0338 tau(i,j) = SQRT(windStress*atmrho)
8053a926ae Jean*0339 #endif
d216b1690c Mart*0340 ENDIF
7311a000e1 Mart*0341
0342
666d41f610 Mart*0343 rhn = (exf_one-stable)*cstanton_1 + stable*cstanton_2
d216b1690c Mart*0344 ren = cDalton
7311a000e1 Mart*0345
d216b1690c Mart*0346 tstar(i,j)=rhn*deltap(i,j)
0347 qstar(i,j)=ren*delq(i,j)
0b902863ee Mart*0348
d216b1690c Mart*0349 ELSE
e1fb02e8f0 Jean*0350
0b902863ee Mart*0351 tstar (i,j) = 0. _d 0
0352 qstar (i,j) = 0. _d 0
0353 ustar (i,j) = 0. _d 0
0354 tau (i,j) = 0. _d 0
0355 rdn (i,j) = 0. _d 0
d216b1690c Mart*0356 ENDIF
0357 ENDDO
0358 ENDDO
0359 DO iter=1,niter_bulk
0360 DO j = 1,sNy
0361 DO i = 1,sNx
0362 IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
0363
0364
bdec91d862 Patr*0365 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0366 ikey_2 = i + (j-1)*sNx + (iter-1)*sNx*sNy
0367 & + (ikey-1)*sNx*sNy*niter_bulk
d216b1690c Mart*0368
0369
0370
0371
0372
0373
bdec91d862 Patr*0374 #endif
0375
d216b1690c Mart*0376 t0 = atemp(i,j,bi,bj)*
0b902863ee Mart*0377 & (exf_one + humid_fac*aqh(i,j,bi,bj))
d216b1690c Mart*0378 huol = ( tstar(i,j)/t0 +
0379 & qstar(i,j)/(exf_one/humid_fac+aqh(i,j,bi,bj))
0380 & )*czol/(ustar(i,j)*ustar(i,j))
0381 #ifdef ALLOW_BULK_LARGEYEAGER04
0382 tmpbulk = MIN(abs(huol),10. _d 0)
0383 huol = SIGN(tmpbulk , huol)
0384 #else
0385
0b902863ee Mart*0386 huol = max(huol,zolmin)
d216b1690c Mart*0387 #endif /* ALLOW_BULK_LARGEYEAGER04 */
0b902863ee Mart*0388 htol = huol*ht/hu
0389 hqol = huol*hq/hu
666d41f610 Mart*0390 stable = exf_half + sign(exf_half, huol)
d216b1690c Mart*0391
8053a926ae Jean*0392
d216b1690c Mart*0393 IF ( solve4Stress ) THEN
0394 #ifdef ALLOW_BULK_LARGEYEAGER04
0395
0396 xsq = SQRT( ABS(exf_one - huol*16. _d 0) )
0397 #else
0398
423768d890 Jean*0399 xsq = MAX(SQRT(ABS(exf_one - 16.*huol)),exf_one)
d216b1690c Mart*0400 #endif /* ALLOW_BULK_LARGEYEAGER04 */
0401 x = SQRT(xsq)
0402 psimh = -psim_fac*huol*stable
0403 & +(exf_one-stable)*
0404 & ( LOG( (exf_one + exf_two*x + xsq)*
0405 & (exf_one+xsq)*.125 _d 0 )
0406 & -exf_two*ATAN(x) + exf_half*pi )
ddea17ebdb Gael*0407 ELSE
0408 psimh = 0. _d 0
d216b1690c Mart*0409 ENDIF
0410 #ifdef ALLOW_BULK_LARGEYEAGER04
0411
0412 xsq = SQRT( ABS(exf_one - htol*16. _d 0) )
0413 #else
0414
423768d890 Jean*0415 xsq = MAX(SQRT(ABS(exf_one - 16.*htol)),exf_one)
d216b1690c Mart*0416 #endif /* ALLOW_BULK_LARGEYEAGER04 */
666d41f610 Mart*0417 psixh = -psim_fac*htol*stable + (exf_one-stable)*
0418 & ( exf_two*LOG(exf_half*(exf_one+xsq)) )
d216b1690c Mart*0419
0420 IF ( solve4Stress ) THEN
7448700841 Mart*0421 #ifdef ALLOW_AUTODIFF_TAMC
8ed910e28e Patr*0422
7448700841 Mart*0423 #endif
d216b1690c Mart*0424
0425 #ifdef ALLOW_BULK_LARGEYEAGER04
0426
0427 dzTmp = (zwln-psimh)/karman
0428 usn = wspeed(i,j,bi,bj)/(exf_one + rdn(i,j)*dzTmp )
0429 #else
0430
666d41f610 Mart*0431
d216b1690c Mart*0432
8053a926ae Jean*0433
0434
d216b1690c Mart*0435 usn = sh(i,j,bi,bj)/(exf_one - rdn(i,j)/karman*psimh)
0436 #endif /* ALLOW_BULK_LARGEYEAGER04 */
0437 usm = MAX(usn, umin)
0438
0439
cc60455fbb Mart*0440 #ifdef ALLOW_DRAG_LARGEYEAGER09
0441
0442 tmpbulk = cdrag_1/usm + cdrag_2 + cdrag_3*usm
0443 & + cdrag_8 * usm**6
0444 tmpbulk = exf_scal_BulkCdn * (
0445 & ( halfRL - SIGN(halfRL, usm-umax) )*tmpbulk
0446 & + ( halfRL + SIGN(halfRL, usm-umax) )*cdragMax
0447 & )
0448 #else
d216b1690c Mart*0449 tmpbulk = exf_scal_BulkCdn *
0450 & ( cdrag_1/usm + cdrag_2 + cdrag_3*usm )
cc60455fbb Mart*0451 #endif
d216b1690c Mart*0452 rdn(i,j) = SQRT(tmpbulk)
0453 #ifdef ALLOW_BULK_LARGEYEAGER04
0454 rd(i,j) = rdn(i,j)/(exf_one + rdn(i,j)*dzTmp)
0455 #else
0456 rd(i,j) = rdn(i,j)/(exf_one - rdn(i,j)/karman*psimh)
0457 #endif /* ALLOW_BULK_LARGEYEAGER04 */
0458 ustar(i,j) = rd(i,j)*sh(i,j,bi,bj)
0459
8053a926ae Jean*0460 #ifdef EXF_CALC_ATMRHO
0461 tau(i,j) = atmrho_loc(i,j)*rd(i,j)*wspeed(i,j,bi,bj)
0462 #else
d216b1690c Mart*0463 tau(i,j) = atmrho*rd(i,j)*wspeed(i,j,bi,bj)
8053a926ae Jean*0464 #endif
d216b1690c Mart*0465 ENDIF
0466
0467
666d41f610 Mart*0468 rhn = (exf_one-stable)*cstanton_1 + stable*cstanton_2
d216b1690c Mart*0469 ren = cDalton
0470
0471
0472 rh = rhn/(exf_one + rhn*(ztln-psixh)/karman)
0473 re = ren/(exf_one + ren*(ztln-psixh)/karman)
bdec91d862 Patr*0474
d216b1690c Mart*0475
0b902863ee Mart*0476 qstar(i,j) = re*delq(i,j)
0477 tstar(i,j) = rh*deltap(i,j)
f1644eaf0d Patr*0478
d216b1690c Mart*0479 ENDIF
0480 ENDDO
0481 ENDDO
8053a926ae Jean*0482
d216b1690c Mart*0483 ENDDO
0484 DO j = 1,sNy
0485 DO i = 1,sNx
0486 IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
bdec91d862 Patr*0487
0488 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0489 ikey_1 = i + (j-1)*sNx + (ikey-1)*sNx*sNy
7311a000e1 Mart*0490
0491
d216b1690c Mart*0492
0493
0494 #endif /* ALLOW_AUTODIFF_TAMC */
0495
0496
666d41f610 Mart*0497 hs(i,j,bi,bj) = atmcp*tau(i,j)*tstar(i,j)
0498 hl(i,j,bi,bj) = flamb*tau(i,j)*qstar(i,j)
bdec91d862 Patr*0499 #ifndef EXF_READ_EVAP
d216b1690c Mart*0500
0501
666d41f610 Mart*0502 evap(i,j,bi,bj) = -recip_rhoConstFresh*tau(i,j)*qstar(i,j)
d216b1690c Mart*0503
0504
bdec91d862 Patr*0505 #endif
56d13a40ed Mart*0506 IF ( useAtmWind .AND. useRelativeWind ) THEN
0507 tmpbulk = tau(i,j)*rd(i,j)
0508 ustress(i,j,bi,bj) = tmpbulk* ( uwind(i,j,bi,bj) -
0509 & 0.5 _d 0 * (uVel(i,j,ks,bi,bj)+uVel(i+1,j,ks,bi,bj)) )
0510 vstress(i,j,bi,bj) = tmpbulk* ( vwind(i,j,bi,bj) -
0511 & 0.5 _d 0 * (vVel(i,j,ks,bi,bj)+vVel(i,j+1,ks,bi,bj)) )
0512 ELSEIF ( useAtmWind ) THEN
423768d890 Jean*0513
0514
d216b1690c Mart*0515
0516
0517
423768d890 Jean*0518 tmpbulk = tau(i,j)*rd(i,j)
0519 ustress(i,j,bi,bj) = tmpbulk*uwind(i,j,bi,bj)
0520 vstress(i,j,bi,bj) = tmpbulk*vwind(i,j,bi,bj)
0521 ENDIF
7311a000e1 Mart*0522
0523
423768d890 Jean*0524 ELSE
0525 IF ( useAtmWind ) ustress(i,j,bi,bj) = 0. _d 0
0526 IF ( useAtmWind ) vstress(i,j,bi,bj) = 0. _d 0
0b902863ee Mart*0527 hflux (i,j,bi,bj) = 0. _d 0
d216b1690c Mart*0528 evap (i,j,bi,bj) = 0. _d 0
0529 hs (i,j,bi,bj) = 0. _d 0
0530 hl (i,j,bi,bj) = 0. _d 0
7311a000e1 Mart*0531
0532
423768d890 Jean*0533 ENDIF
0534
d216b1690c Mart*0535 ENDDO
0536 ENDDO
0537 ENDDO
0538 ENDDO
bdec91d862 Patr*0539
423768d890 Jean*0540 #endif /* ALLOW_ATM_TEMP */
7109a141b2 Patr*0541 #endif /* ALLOW_BULKFORMULAE */
0542
d216b1690c Mart*0543 RETURN
0544 END