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