| 1 | molod | 1.2 | C $Header: /u/gcmpack/MITgcm_contrib/plumes/external_forcing.F,v 1.1 2004/05/13 22:21:45 molod Exp $ | 
| 2 | molod | 1.1 | C $Name:  $ | 
| 3 |  |  |  | 
| 4 |  |  | #include "PACKAGES_CONFIG.h" | 
| 5 |  |  | #include "CPP_OPTIONS.h" | 
| 6 |  |  | #ifdef ALLOW_OBCS | 
| 7 |  |  | # include "OBCS_OPTIONS.h" | 
| 8 |  |  | #endif | 
| 9 |  |  |  | 
| 10 |  |  | CBOP | 
| 11 |  |  | C     !ROUTINE: EXTERNAL_FORCING_U | 
| 12 |  |  | C     !INTERFACE: | 
| 13 |  |  | SUBROUTINE EXTERNAL_FORCING_U( | 
| 14 |  |  | I           iMin, iMax, jMin, jMax,bi,bj,kLev, | 
| 15 |  |  | I           myCurrentTime,myThid) | 
| 16 |  |  | C     !DESCRIPTION: \bv | 
| 17 |  |  | C     *==========================================================* | 
| 18 |  |  | C     | S/R EXTERNAL_FORCING_U | 
| 19 |  |  | C     | o Contains problem specific forcing for zonal velocity. | 
| 20 |  |  | C     *==========================================================* | 
| 21 |  |  | C     | Adds terms to gU for forcing by external sources | 
| 22 |  |  | C     | e.g. wind stress, bottom friction etc.................. | 
| 23 |  |  | C     *==========================================================* | 
| 24 |  |  | C     \ev | 
| 25 |  |  |  | 
| 26 |  |  | C     !USES: | 
| 27 |  |  | IMPLICIT NONE | 
| 28 |  |  | C     == Global data == | 
| 29 |  |  | #include "SIZE.h" | 
| 30 |  |  | #include "EEPARAMS.h" | 
| 31 |  |  | #include "PARAMS.h" | 
| 32 |  |  | #include "GRID.h" | 
| 33 |  |  | #include "DYNVARS.h" | 
| 34 |  |  | #include "FFIELDS.h" | 
| 35 |  |  |  | 
| 36 |  |  | C     !INPUT/OUTPUT PARAMETERS: | 
| 37 |  |  | C     == Routine arguments == | 
| 38 |  |  | C     iMin - Working range of tile for applying forcing. | 
| 39 |  |  | C     iMax | 
| 40 |  |  | C     jMin | 
| 41 |  |  | C     jMax | 
| 42 |  |  | C     kLev | 
| 43 |  |  | INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj | 
| 44 |  |  | _RL myCurrentTime | 
| 45 |  |  | INTEGER myThid | 
| 46 |  |  |  | 
| 47 |  |  | C     !LOCAL VARIABLES: | 
| 48 |  |  | C     == Local variables == | 
| 49 |  |  | C     Loop counters | 
| 50 |  |  | INTEGER I, J | 
| 51 |  |  | C     number of surface interface layer | 
| 52 |  |  | INTEGER kSurface | 
| 53 |  |  | CEOP | 
| 54 |  |  |  | 
| 55 |  |  | if ( buoyancyRelation .eq. 'ATMOSPHERIC' ) then | 
| 56 |  |  | kSurface = 0 | 
| 57 |  |  | elseif ( buoyancyRelation .eq. 'OCEANICP' ) then | 
| 58 |  |  | kSurface = Nr | 
| 59 |  |  | else | 
| 60 |  |  | kSurface = 1 | 
| 61 |  |  | endif | 
| 62 |  |  |  | 
| 63 |  |  | C--   Forcing term | 
| 64 |  |  | #ifdef ALLOW_AIM | 
| 65 |  |  | IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U( | 
| 66 |  |  | &                      iMin,iMax, jMin,jMax, bi,bj, kLev, | 
| 67 |  |  | &                      myCurrentTime, myThid ) | 
| 68 |  |  | #endif /* ALLOW_AIM */ | 
| 69 |  |  | C AMM | 
| 70 |  |  | #ifdef ALLOW_FIZHI | 
| 71 |  |  | IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U( | 
| 72 |  |  | &                      iMin,iMax, jMin,jMax, bi,bj, kLev, | 
| 73 |  |  | &                      myCurrentTime, myThid ) | 
| 74 |  |  | #endif /* ALLOW_FIZHI */ | 
| 75 |  |  | C AMM | 
| 76 |  |  |  | 
| 77 | molod | 1.2 | C AMM | 
| 78 |  |  | CALL PLUME_TENDENCY_APPLY_U( | 
| 79 |  |  | &                      iMin,iMax, jMin,jMax, bi,bj, kLev, | 
| 80 |  |  | &                      myCurrentTime, myThid ) | 
| 81 |  |  | C AMM | 
| 82 |  |  |  | 
| 83 | molod | 1.1 | C     Add windstress momentum impulse into the top-layer | 
| 84 |  |  | IF ( kLev .EQ. kSurface ) THEN | 
| 85 |  |  | DO j=jMin,jMax | 
| 86 |  |  | DO i=iMin,iMax | 
| 87 |  |  | gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj) | 
| 88 |  |  | &   +foFacMom*surfaceTendencyU(i,j,bi,bj) | 
| 89 |  |  | &   *_maskW(i,j,kLev,bi,bj) | 
| 90 |  |  | ENDDO | 
| 91 |  |  | ENDDO | 
| 92 |  |  | ENDIF | 
| 93 |  |  |  | 
| 94 |  |  | #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE)) | 
| 95 |  |  | IF (useOBCS) THEN | 
| 96 |  |  | CALL OBCS_SPONGE_U( | 
| 97 |  |  | I           iMin, iMax, jMin, jMax,bi,bj,kLev, | 
| 98 |  |  | I           myCurrentTime,myThid) | 
| 99 |  |  | ENDIF | 
| 100 |  |  | #endif | 
| 101 |  |  |  | 
| 102 |  |  | RETURN | 
| 103 |  |  | END | 
| 104 |  |  | CBOP | 
| 105 |  |  | C     !ROUTINE: EXTERNAL_FORCING_V | 
| 106 |  |  | C     !INTERFACE: | 
| 107 |  |  | SUBROUTINE EXTERNAL_FORCING_V( | 
| 108 |  |  | I           iMin, iMax, jMin, jMax,bi,bj,kLev, | 
| 109 |  |  | I           myCurrentTime,myThid) | 
| 110 |  |  | C     !DESCRIPTION: \bv | 
| 111 |  |  | C     *==========================================================* | 
| 112 |  |  | C     | S/R EXTERNAL_FORCING_V | 
| 113 |  |  | C     | o Contains problem specific forcing for merid velocity. | 
| 114 |  |  | C     *==========================================================* | 
| 115 |  |  | C     | Adds terms to gV for forcing by external sources | 
| 116 |  |  | C     | e.g. wind stress, bottom friction etc.................. | 
| 117 |  |  | C     *==========================================================* | 
| 118 |  |  | C     \ev | 
| 119 |  |  |  | 
| 120 |  |  | C     !USES: | 
| 121 |  |  | IMPLICIT NONE | 
| 122 |  |  | C     == Global data == | 
| 123 |  |  | #include "SIZE.h" | 
| 124 |  |  | #include "EEPARAMS.h" | 
| 125 |  |  | #include "PARAMS.h" | 
| 126 |  |  | #include "GRID.h" | 
| 127 |  |  | #include "DYNVARS.h" | 
| 128 |  |  | #include "FFIELDS.h" | 
| 129 |  |  |  | 
| 130 |  |  | C     !INPUT/OUTPUT PARAMETERS: | 
| 131 |  |  | C     == Routine arguments == | 
| 132 |  |  | C     iMin - Working range of tile for applying forcing. | 
| 133 |  |  | C     iMax | 
| 134 |  |  | C     jMin | 
| 135 |  |  | C     jMax | 
| 136 |  |  | C     kLev | 
| 137 |  |  | INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj | 
| 138 |  |  | _RL myCurrentTime | 
| 139 |  |  | INTEGER myThid | 
| 140 |  |  |  | 
| 141 |  |  | C     !LOCAL VARIABLES: | 
| 142 |  |  | C     == Local variables == | 
| 143 |  |  | C     Loop counters | 
| 144 |  |  | INTEGER I, J | 
| 145 |  |  | C     number of surface interface layer | 
| 146 |  |  | INTEGER kSurface | 
| 147 |  |  | CEOP | 
| 148 |  |  |  | 
| 149 |  |  | if ( buoyancyRelation .eq. 'ATMOSPHERIC' ) then | 
| 150 |  |  | kSurface = 0 | 
| 151 |  |  | elseif ( buoyancyRelation .eq. 'OCEANICP' ) then | 
| 152 |  |  | kSurface = Nr | 
| 153 |  |  | else | 
| 154 |  |  | kSurface = 1 | 
| 155 |  |  | endif | 
| 156 |  |  |  | 
| 157 |  |  | C--   Forcing term | 
| 158 |  |  | #ifdef ALLOW_AIM | 
| 159 |  |  | IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V( | 
| 160 |  |  | &                      iMin,iMax, jMin,jMax, bi,bj, kLev, | 
| 161 |  |  | &                      myCurrentTime, myThid ) | 
| 162 |  |  | #endif /* ALLOW_AIM */ | 
| 163 |  |  |  | 
| 164 |  |  | C AMM | 
| 165 |  |  | #ifdef ALLOW_FIZHI | 
| 166 |  |  | IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V( | 
| 167 |  |  | &                      iMin,iMax, jMin,jMax, bi,bj, kLev, | 
| 168 |  |  | &                      myCurrentTime, myThid ) | 
| 169 |  |  | #endif /* ALLOW_FIZHI */ | 
| 170 |  |  | C AMM | 
| 171 | molod | 1.2 |  | 
| 172 |  |  | C AMM | 
| 173 |  |  | CALL PLUME_TENDENCY_APPLY_V( | 
| 174 |  |  | &                      iMin,iMax, jMin,jMax, bi,bj, kLev, | 
| 175 |  |  | &                      myCurrentTime, myThid ) | 
| 176 |  |  | C AMM | 
| 177 |  |  |  | 
| 178 | molod | 1.1 | C     Add windstress momentum impulse into the top-layer | 
| 179 |  |  | IF ( kLev .EQ. kSurface ) THEN | 
| 180 |  |  | DO j=jMin,jMax | 
| 181 |  |  | DO i=iMin,iMax | 
| 182 |  |  | gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj) | 
| 183 |  |  | &   +foFacMom*surfaceTendencyV(i,j,bi,bj) | 
| 184 |  |  | &   *_maskS(i,j,kLev,bi,bj) | 
| 185 |  |  | ENDDO | 
| 186 |  |  | ENDDO | 
| 187 |  |  | ENDIF | 
| 188 |  |  |  | 
| 189 |  |  | #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE)) | 
| 190 |  |  | IF (useOBCS) THEN | 
| 191 |  |  | CALL OBCS_SPONGE_V( | 
| 192 |  |  | I           iMin, iMax, jMin, jMax,bi,bj,kLev, | 
| 193 |  |  | I           myCurrentTime,myThid) | 
| 194 |  |  | ENDIF | 
| 195 |  |  | #endif | 
| 196 |  |  |  | 
| 197 |  |  | RETURN | 
| 198 |  |  | END | 
| 199 |  |  | CBOP | 
| 200 |  |  | C     !ROUTINE: EXTERNAL_FORCING_T | 
| 201 |  |  | C     !INTERFACE: | 
| 202 |  |  | SUBROUTINE EXTERNAL_FORCING_T( | 
| 203 |  |  | I           iMin, iMax, jMin, jMax,bi,bj,kLev, | 
| 204 |  |  | I           myCurrentTime,myThid) | 
| 205 |  |  | C     !DESCRIPTION: \bv | 
| 206 |  |  | C     *==========================================================* | 
| 207 |  |  | C     | S/R EXTERNAL_FORCING_T | 
| 208 |  |  | C     | o Contains problem specific forcing for temperature. | 
| 209 |  |  | C     *==========================================================* | 
| 210 |  |  | C     | Adds terms to gT for forcing by external sources | 
| 211 |  |  | C     | e.g. heat flux, climatalogical relaxation.............. | 
| 212 |  |  | C     *==========================================================* | 
| 213 |  |  | C     \ev | 
| 214 |  |  |  | 
| 215 |  |  | C     !USES: | 
| 216 |  |  | IMPLICIT NONE | 
| 217 |  |  | C     == Global data == | 
| 218 |  |  | #include "SIZE.h" | 
| 219 |  |  | #include "EEPARAMS.h" | 
| 220 |  |  | #include "PARAMS.h" | 
| 221 |  |  | #include "GRID.h" | 
| 222 |  |  | #include "DYNVARS.h" | 
| 223 |  |  | #include "FFIELDS.h" | 
| 224 |  |  |  | 
| 225 |  |  | C     !INPUT/OUTPUT PARAMETERS: | 
| 226 |  |  | C     == Routine arguments == | 
| 227 |  |  | C     iMin - Working range of tile for applying forcing. | 
| 228 |  |  | C     iMax | 
| 229 |  |  | C     jMin | 
| 230 |  |  | C     jMax | 
| 231 |  |  | C     kLev | 
| 232 |  |  | INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj | 
| 233 |  |  | _RL myCurrentTime | 
| 234 |  |  | INTEGER myThid | 
| 235 |  |  | CEndOfInterface | 
| 236 |  |  |  | 
| 237 |  |  | C     !LOCAL VARIABLES: | 
| 238 |  |  | C     == Local variables == | 
| 239 |  |  | C     Loop counters | 
| 240 |  |  | INTEGER I, J | 
| 241 |  |  | C     number of surface interface layer | 
| 242 |  |  | INTEGER kSurface | 
| 243 |  |  | #ifdef SHORTWAVE_HEATING | 
| 244 |  |  | integer two | 
| 245 |  |  | _RL minusone | 
| 246 |  |  | parameter (two=2,minusone=-1.) | 
| 247 |  |  | _RL swfracb(two) | 
| 248 |  |  | INTEGER kp1 | 
| 249 |  |  | #endif | 
| 250 |  |  | CEOP | 
| 251 |  |  |  | 
| 252 |  |  | if ( buoyancyRelation .eq. 'ATMOSPHERIC' ) then | 
| 253 |  |  | kSurface = 0 | 
| 254 |  |  | elseif ( buoyancyRelation .eq. 'OCEANICP' ) then | 
| 255 |  |  | kSurface = Nr | 
| 256 |  |  | else | 
| 257 |  |  | kSurface = 1 | 
| 258 |  |  | endif | 
| 259 |  |  |  | 
| 260 |  |  | C--   Forcing term | 
| 261 |  |  | #ifdef ALLOW_AIM | 
| 262 |  |  | IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T( | 
| 263 |  |  | &                      iMin,iMax, jMin,jMax, bi,bj, kLev, | 
| 264 |  |  | &                      myCurrentTime, myThid ) | 
| 265 |  |  | #endif /* ALLOW_AIM */ | 
| 266 |  |  |  | 
| 267 |  |  | C AMM | 
| 268 |  |  | #ifdef ALLOW_FIZHI | 
| 269 |  |  | IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T( | 
| 270 |  |  | &                      iMin,iMax, jMin,jMax, bi,bj, kLev, | 
| 271 |  |  | &                      myCurrentTime, myThid ) | 
| 272 |  |  | #endif /* ALLOW_FIZHI */ | 
| 273 |  |  | C AMM | 
| 274 | molod | 1.2 |  | 
| 275 |  |  | C AMM | 
| 276 |  |  | CALL PLUME_TENDENCY_APPLY_T( | 
| 277 |  |  | &                      iMin,iMax, jMin,jMax, bi,bj, kLev, | 
| 278 |  |  | &                      myCurrentTime, myThid ) | 
| 279 |  |  | C AMM | 
| 280 |  |  |  | 
| 281 | molod | 1.1 | C     Add heat in top-layer | 
| 282 |  |  | IF ( kLev .EQ. kSurface ) THEN | 
| 283 |  |  | DO j=jMin,jMax | 
| 284 |  |  | DO i=iMin,iMax | 
| 285 |  |  | gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj) | 
| 286 |  |  | &     +maskC(i,j,kLev,bi,bj)*surfaceTendencyT(i,j,bi,bj) | 
| 287 |  |  | ENDDO | 
| 288 |  |  | ENDDO | 
| 289 |  |  | ENDIF | 
| 290 |  |  |  | 
| 291 |  |  | #ifdef SHORTWAVE_HEATING | 
| 292 |  |  | C Penetrating SW radiation | 
| 293 |  |  | kp1 = klev+1 | 
| 294 |  |  | swfracb(1)=abs(rF(klev)) | 
| 295 |  |  | swfracb(2)=abs(rF(klev+1)) | 
| 296 |  |  | CALL SWFRAC( | 
| 297 |  |  | I     two,minusone, | 
| 298 |  |  | I     myCurrentTime,myThid, | 
| 299 |  |  | U     swfracb) | 
| 300 |  |  | IF (klev.EQ.Nr) THEN | 
| 301 |  |  | kp1 = klev | 
| 302 |  |  | swfracb(2)=0. _d 0 | 
| 303 |  |  | ENDIF | 
| 304 |  |  | DO j=jMin,jMax | 
| 305 |  |  | DO i=iMin,iMax | 
| 306 |  |  | gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj) | 
| 307 |  |  | &   -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj) | 
| 308 |  |  | &                   -swfracb(2)*maskC(i,j,kp1, bi,bj)) | 
| 309 |  |  | &    *recip_Cp*recip_rhoConst*recip_drF(klev) | 
| 310 |  |  | ENDDO | 
| 311 |  |  | ENDDO | 
| 312 |  |  | #endif | 
| 313 |  |  |  | 
| 314 |  |  | #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE)) | 
| 315 |  |  | IF (useOBCS) THEN | 
| 316 |  |  | CALL OBCS_SPONGE_T( | 
| 317 |  |  | I           iMin, iMax, jMin, jMax,bi,bj,kLev, | 
| 318 |  |  | I           myCurrentTime,myThid) | 
| 319 |  |  | ENDIF | 
| 320 |  |  | #endif | 
| 321 |  |  |  | 
| 322 |  |  | RETURN | 
| 323 |  |  | END | 
| 324 |  |  | CBOP | 
| 325 |  |  | C     !ROUTINE: EXTERNAL_FORCING_S | 
| 326 |  |  | C     !INTERFACE: | 
| 327 |  |  | SUBROUTINE EXTERNAL_FORCING_S( | 
| 328 |  |  | I           iMin, iMax, jMin, jMax,bi,bj,kLev, | 
| 329 |  |  | I           myCurrentTime,myThid) | 
| 330 |  |  |  | 
| 331 |  |  | C     !DESCRIPTION: \bv | 
| 332 |  |  | C     *==========================================================* | 
| 333 |  |  | C     | S/R EXTERNAL_FORCING_S | 
| 334 |  |  | C     | o Contains problem specific forcing for merid velocity. | 
| 335 |  |  | C     *==========================================================* | 
| 336 |  |  | C     | Adds terms to gS for forcing by external sources | 
| 337 |  |  | C     | e.g. fresh-water flux, climatalogical relaxation....... | 
| 338 |  |  | C     *==========================================================* | 
| 339 |  |  | C     \ev | 
| 340 |  |  |  | 
| 341 |  |  | C     !USES: | 
| 342 |  |  | IMPLICIT NONE | 
| 343 |  |  | C     == Global data == | 
| 344 |  |  | #include "SIZE.h" | 
| 345 |  |  | #include "EEPARAMS.h" | 
| 346 |  |  | #include "PARAMS.h" | 
| 347 |  |  | #include "GRID.h" | 
| 348 |  |  | #include "DYNVARS.h" | 
| 349 |  |  | #include "FFIELDS.h" | 
| 350 |  |  |  | 
| 351 |  |  | C     !INPUT/OUTPUT PARAMETERS: | 
| 352 |  |  | C     == Routine arguments == | 
| 353 |  |  | C     iMin - Working range of tile for applying forcing. | 
| 354 |  |  | C     iMax | 
| 355 |  |  | C     jMin | 
| 356 |  |  | C     jMax | 
| 357 |  |  | C     kLev | 
| 358 |  |  | INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj | 
| 359 |  |  | _RL myCurrentTime | 
| 360 |  |  | INTEGER myThid | 
| 361 |  |  |  | 
| 362 |  |  | C     !LOCAL VARIABLES: | 
| 363 |  |  | C     == Local variables == | 
| 364 |  |  | C     Loop counters | 
| 365 |  |  | INTEGER I, J | 
| 366 |  |  | C     number of surface interface layer | 
| 367 |  |  | INTEGER kSurface | 
| 368 |  |  | CEOP | 
| 369 |  |  |  | 
| 370 |  |  | if ( buoyancyRelation .eq. 'ATMOSPHERIC' ) then | 
| 371 |  |  | kSurface = 0 | 
| 372 |  |  | elseif ( buoyancyRelation .eq. 'OCEANICP' ) then | 
| 373 |  |  | kSurface = Nr | 
| 374 |  |  | else | 
| 375 |  |  | kSurface = 1 | 
| 376 |  |  | endif | 
| 377 |  |  |  | 
| 378 |  |  |  | 
| 379 |  |  | C--   Forcing term | 
| 380 |  |  | #ifdef ALLOW_AIM | 
| 381 |  |  | IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S( | 
| 382 |  |  | &                      iMin,iMax, jMin,jMax, bi,bj, kLev, | 
| 383 |  |  | &                      myCurrentTime, myThid ) | 
| 384 |  |  | #endif /* ALLOW_AIM */ | 
| 385 |  |  |  | 
| 386 |  |  | C AMM | 
| 387 |  |  | #ifdef ALLOW_FIZHI | 
| 388 |  |  | IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S( | 
| 389 |  |  | &                      iMin,iMax, jMin,jMax, bi,bj, kLev, | 
| 390 |  |  | &                      myCurrentTime, myThid ) | 
| 391 |  |  | #endif /* ALLOW_FIZHI */ | 
| 392 |  |  | C AMM | 
| 393 |  |  | C     Add fresh-water in top-layer | 
| 394 |  |  | IF ( kLev .EQ. kSurface ) THEN | 
| 395 |  |  | DO j=jMin,jMax | 
| 396 |  |  | DO i=iMin,iMax | 
| 397 |  |  | gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj) | 
| 398 |  |  | &   +maskC(i,j,kLev,bi,bj)*surfaceTendencyS(i,j,bi,bj) | 
| 399 |  |  | ENDDO | 
| 400 |  |  | ENDDO | 
| 401 |  |  | ENDIF | 
| 402 |  |  |  | 
| 403 |  |  | #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE)) | 
| 404 |  |  | IF (useOBCS) THEN | 
| 405 |  |  | CALL OBCS_SPONGE_S( | 
| 406 |  |  | I           iMin, iMax, jMin, jMax,bi,bj,kLev, | 
| 407 |  |  | I           myCurrentTime,myThid) | 
| 408 |  |  | ENDIF | 
| 409 |  |  | #endif | 
| 410 |  |  |  | 
| 411 |  |  | RETURN | 
| 412 |  |  | END |