[3100] | 1 | ! |
---|
| 2 | MODULE slab_heat_transp_mod |
---|
| 3 | ! |
---|
| 4 | ! Slab ocean : temperature tendencies due to horizontal diffusion |
---|
| 5 | ! and / or Ekman transport |
---|
| 6 | |
---|
| 7 | USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo |
---|
| 8 | IMPLICIT NONE |
---|
| 9 | |
---|
| 10 | ! Variables copied over from dyn3d dynamics: |
---|
| 11 | REAL,SAVE,ALLOCATABLE :: fext(:) ! Coriolis f times cell area |
---|
| 12 | !$OMP THREADPRIVATE(fext) |
---|
| 13 | REAL,SAVE,ALLOCATABLE :: beta(:) ! df/dy |
---|
| 14 | !$OMP THREADPRIVATE(beta) |
---|
| 15 | REAL,SAVE,ALLOCATABLE :: unsairez(:) ! 1/cell area |
---|
| 16 | !$OMP THREADPRIVATE(unsairez) |
---|
| 17 | REAL,SAVE,ALLOCATABLE :: unsaire(:) |
---|
| 18 | !$OMP THREADPRIVATE(unsaire) |
---|
| 19 | REAL,SAVE,ALLOCATABLE :: cu(:) ! cell longitude dim (m) |
---|
| 20 | !$OMP THREADPRIVATE(cu) |
---|
| 21 | REAL,SAVE,ALLOCATABLE :: cv(:) ! cell latitude dim (m) |
---|
| 22 | !$OMP THREADPRIVATE(cv) |
---|
| 23 | REAL,SAVE,ALLOCATABLE :: cuvsurcv(:) ! cu/cv (v points) |
---|
| 24 | !$OMP THREADPRIVATE(cuvsurcv) |
---|
| 25 | REAL,SAVE,ALLOCATABLE :: cvusurcu(:) ! cv/cu (u points) |
---|
| 26 | !$OMP THREADPRIVATE(cvusurcu) |
---|
| 27 | REAL,SAVE,ALLOCATABLE :: aire(:) ! cell area |
---|
| 28 | !$OMP THREADPRIVATE(aire) |
---|
| 29 | REAL,SAVE :: apoln ! area of north pole points |
---|
| 30 | !$OMP THREADPRIVATE(apoln) |
---|
| 31 | REAL,SAVE :: apols ! area of south pole points |
---|
| 32 | !$OMP THREADPRIVATE(apols) |
---|
| 33 | REAL,SAVE,ALLOCATABLE :: aireu(:) ! area of u cells |
---|
| 34 | !$OMP THREADPRIVATE(aireu) |
---|
| 35 | REAL,SAVE,ALLOCATABLE :: airev(:) ! area of v cells |
---|
| 36 | !$OMP THREADPRIVATE(airev) |
---|
| 37 | |
---|
| 38 | ! Local parameters for slab transport |
---|
| 39 | !LOGICAL,SAVE :: alpha_var ! variable coef for deep temp (1 layer) |
---|
| 40 | !!$OMP THREADPRIVATE(alpha_var) |
---|
| 41 | LOGICAL,SAVE :: slab_upstream ! upstream scheme ? (1 layer) |
---|
| 42 | !$OMP THREADPRIVATE(slab_upstream) |
---|
| 43 | LOGICAL,SAVE :: slab_sverdrup ! use wind stress curl at equator |
---|
| 44 | !$OMP THREADPRIVATE(slab_sverdrup) |
---|
| 45 | LOGICAL,SAVE :: slab_tyeq ! use merid wind stress at equator |
---|
| 46 | !$OMP THREADPRIVATE(slab_tyeq) |
---|
| 47 | LOGICAL,SAVE :: ekman_zonadv ! use zonal advection by Ekman currents |
---|
| 48 | !$OMP THREADPRIVATE(ekman_zonadv) |
---|
| 49 | LOGICAL,SAVE :: ekman_zonavg ! zonally averaged wind stress |
---|
| 50 | !$OMP THREADPRIVATE(ekman_zonavg) |
---|
| 51 | |
---|
| 52 | REAL,SAVE :: alpham |
---|
| 53 | !$OMP THREADPRIVATE(alpham) |
---|
| 54 | REAL,SAVE :: gmkappa ! GM scheme coefficient |
---|
| 55 | !$OMP THREADPRIVATE(gmkappa) |
---|
| 56 | REAL,SAVE :: gm_smax ! GM scheme limiter (max slope) |
---|
| 57 | !$OMP THREADPRIVATE(gm_smax) |
---|
| 58 | |
---|
| 59 | ! geometry variables : f, beta, mask... |
---|
| 60 | REAL,SAVE,ALLOCATABLE :: zmasqu(:) ! continent mask for zonal mass flux |
---|
| 61 | !$OMP THREADPRIVATE(zmasqu) |
---|
| 62 | REAL,SAVE,ALLOCATABLE :: zmasqv(:) ! continent mask for merid mass flux |
---|
| 63 | !$OMP THREADPRIVATE(zmasqv) |
---|
| 64 | REAL,SAVE,ALLOCATABLE :: unsfv(:) ! 1/f, v points |
---|
| 65 | !$OMP THREADPRIVATE(unsfv) |
---|
| 66 | REAL,SAVE,ALLOCATABLE :: unsbv(:) ! 1/beta |
---|
| 67 | !$OMP THREADPRIVATE(unsbv) |
---|
| 68 | REAL,SAVE,ALLOCATABLE :: unsev(:) ! 1/epsilon (drag) |
---|
| 69 | !$OMP THREADPRIVATE(unsev) |
---|
| 70 | REAL,SAVE,ALLOCATABLE :: unsfu(:) ! 1/F, u points |
---|
| 71 | !$OMP THREADPRIVATE(unsfu) |
---|
| 72 | REAL,SAVE,ALLOCATABLE :: unseu(:) |
---|
| 73 | !$OMP THREADPRIVATE(unseu) |
---|
| 74 | |
---|
| 75 | ! Routines from dyn3d, valid on global dynamics grid only: |
---|
| 76 | PRIVATE :: gr_fi_dyn, gr_dyn_fi ! to go between 1D nd 2D horiz grid |
---|
| 77 | PRIVATE :: gr_scal_v,gr_v_scal,gr_scal_u ! change on 2D grid U,V, T points |
---|
| 78 | PRIVATE :: grad,diverg |
---|
| 79 | |
---|
| 80 | CONTAINS |
---|
| 81 | |
---|
| 82 | SUBROUTINE ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez_,fext_,unsaire_,& |
---|
| 83 | cu_,cuvsurcv_,cv_,cvusurcu_, & |
---|
| 84 | aire_,apoln_,apols_, & |
---|
| 85 | aireu_,airev_,rlatv) |
---|
| 86 | ! USE comcstfi_mod, ONLY: omeg, rad |
---|
| 87 | ! number of points in lon, lat |
---|
| 88 | IMPLICIT NONE |
---|
| 89 | ! Routine copies some geometry variables from the dynamical core |
---|
| 90 | ! see global vars for meaning |
---|
| 91 | INTEGER,INTENT(IN) :: ip1jm |
---|
| 92 | INTEGER,INTENT(IN) :: ip1jmp1 |
---|
| 93 | REAL,INTENT(IN) :: unsairez_(ip1jm) |
---|
| 94 | REAL,INTENT(IN) :: fext_(ip1jm) |
---|
| 95 | REAL,INTENT(IN) :: unsaire_(ip1jmp1) |
---|
| 96 | REAL,INTENT(IN) :: cu_(ip1jmp1) |
---|
| 97 | REAL,INTENT(IN) :: cuvsurcv_(ip1jm) |
---|
| 98 | REAL,INTENT(IN) :: cv_(ip1jm) |
---|
| 99 | REAL,INTENT(IN) :: cvusurcu_(ip1jmp1) |
---|
| 100 | REAL,INTENT(IN) :: aire_(ip1jmp1) |
---|
| 101 | REAL,INTENT(IN) :: apoln_ |
---|
| 102 | REAL,INTENT(IN) :: apols_ |
---|
| 103 | REAL,INTENT(IN) :: aireu_(ip1jmp1) |
---|
| 104 | REAL,INTENT(IN) :: airev_(ip1jm) |
---|
| 105 | REAL,INTENT(IN) :: rlatv(nbp_lat-1) |
---|
| 106 | |
---|
| 107 | ! Sanity check on dimensions |
---|
| 108 | if ((ip1jm.ne.((nbp_lon+1)*(nbp_lat-1))).or. & |
---|
| 109 | (ip1jmp1.ne.((nbp_lon+1)*nbp_lat))) then |
---|
| 110 | write(*,*) "ini_slab_transp_geom Error: wrong array sizes" |
---|
| 111 | stop |
---|
| 112 | endif |
---|
| 113 | ! Allocations could be done only on master process/thread... |
---|
| 114 | allocate(unsairez(ip1jm)) |
---|
| 115 | unsairez(:)=unsairez_(:) |
---|
| 116 | allocate(fext(ip1jm)) |
---|
| 117 | fext(:)=fext_(:) |
---|
| 118 | allocate(unsaire(ip1jmp1)) |
---|
| 119 | unsaire(:)=unsaire_(:) |
---|
| 120 | allocate(cu(ip1jmp1)) |
---|
| 121 | cu(:)=cu_(:) |
---|
| 122 | allocate(cuvsurcv(ip1jm)) |
---|
| 123 | cuvsurcv(:)=cuvsurcv_(:) |
---|
| 124 | allocate(cv(ip1jm)) |
---|
| 125 | cv(:)=cv_(:) |
---|
| 126 | allocate(cvusurcu(ip1jmp1)) |
---|
| 127 | cvusurcu(:)=cvusurcu_(:) |
---|
| 128 | allocate(aire(ip1jmp1)) |
---|
| 129 | aire(:)=aire_(:) |
---|
| 130 | apoln=apoln_ |
---|
| 131 | apols=apols_ |
---|
| 132 | allocate(aireu(ip1jmp1)) |
---|
| 133 | aireu(:)=aireu_(:) |
---|
| 134 | allocate(airev(ip1jm)) |
---|
| 135 | airev(:)=airev_(:) |
---|
| 136 | allocate(beta(nbp_lat-1)) |
---|
| 137 | ! beta(:)=2*omeg*cos(rlatv(:))/rad |
---|
| 138 | beta(:)=rlatv(:) ! Temporarily storing rlatv in beta, will be updated in ini_slab_transp |
---|
| 139 | END SUBROUTINE ini_slab_transp_geom |
---|
| 140 | |
---|
| 141 | SUBROUTINE ini_slab_transp(zmasq) |
---|
| 142 | |
---|
| 143 | ! USE ioipsl_getin_p_mod, only: getin_p |
---|
| 144 | USE IOIPSL, ONLY : getin |
---|
| 145 | USE comcstfi_mod, ONLY: omeg, rad |
---|
| 146 | IMPLICIT NONE |
---|
| 147 | |
---|
| 148 | REAL, INTENT (IN) :: zmasq(klon_glo) ! ocean / continent mask, 1=continent |
---|
| 149 | REAL zmasq_2d((nbp_lon+1)*nbp_lat) |
---|
| 150 | REAL ff((nbp_lon+1)*(nbp_lat-1)) ! Coriolis parameter |
---|
| 151 | REAL eps ! epsilon friction timescale (s-1) |
---|
| 152 | ! INTEGER :: slab_ekman |
---|
| 153 | LOGICAL :: slab_ekman |
---|
| 154 | INTEGER i |
---|
| 155 | INTEGER :: iim,iip1,jjp1,ip1jm,ip1jmp1 |
---|
| 156 | |
---|
| 157 | ! Some definition for grid size |
---|
| 158 | ip1jm=(nbp_lon+1)*(nbp_lat-1) |
---|
| 159 | ip1jmp1=(nbp_lon+1)*nbp_lat |
---|
| 160 | iim=nbp_lon |
---|
| 161 | iip1=nbp_lon+1 |
---|
| 162 | jjp1=nbp_lat |
---|
| 163 | ip1jm=(nbp_lon+1)*(nbp_lat-1) |
---|
| 164 | ip1jmp1=(nbp_lon+1)*nbp_lat |
---|
| 165 | |
---|
| 166 | ! Options for Heat transport |
---|
| 167 | ! Alpha variable? (for 1.5-layer Ekman) |
---|
| 168 | ! alpha_var=.FALSE. |
---|
| 169 | ! CALL getin('slab_alphav',alpha_var) |
---|
| 170 | ! print *,'alpha variable',alpha_var |
---|
| 171 | ! centered ou upstream scheme for meridional transport |
---|
| 172 | slab_upstream=.FALSE. |
---|
| 173 | CALL getin('slab_upstream',slab_upstream) |
---|
| 174 | print *,'upstream slab scheme',slab_upstream |
---|
| 175 | ! Sverdrup balance (wind stress curl) at equator ? |
---|
| 176 | slab_sverdrup=.TRUE. |
---|
| 177 | CALL getin('slab_sverdrup',slab_sverdrup) |
---|
| 178 | print *,'Sverdrup balance',slab_sverdrup |
---|
| 179 | ! Use tauy for meridional flux at equator ? |
---|
| 180 | slab_tyeq=.TRUE. |
---|
| 181 | CALL getin('slab_tyeq',slab_tyeq) |
---|
| 182 | print *,'Tauy forcing at equator',slab_tyeq |
---|
| 183 | ! Zonal Ekman transport ? (Or meridional only) |
---|
| 184 | ekman_zonadv=.TRUE. |
---|
| 185 | CALL getin('slab_ekman_zonadv',ekman_zonadv) |
---|
| 186 | print *,'Use Ekman flow in zonal direction',ekman_zonadv |
---|
| 187 | ! Zonally average the Ekman flux ? |
---|
| 188 | ekman_zonavg=.FALSE. |
---|
| 189 | CALL getin('slab_ekman_zonavg',ekman_zonavg) |
---|
| 190 | print *,'Use zonally-averaged wind stress ?',ekman_zonavg |
---|
| 191 | ! Value of alpha (Ekman 1.5-layer Tdeep computation) |
---|
| 192 | alpham=2./3. |
---|
| 193 | CALL getin('slab_alpha',alpham) |
---|
| 194 | print *,'slab_alpha',alpham |
---|
| 195 | ! GM k coefficient (m2/s) for 2-layers |
---|
| 196 | gmkappa=1000. |
---|
| 197 | CALL getin('slab_gmkappa',gmkappa) |
---|
| 198 | print *,'slab_gmkappa',gmkappa |
---|
| 199 | ! GM max slope for 2-layers |
---|
| 200 | gm_smax=2e-3 |
---|
| 201 | CALL getin('slab_gm_smax',gm_smax) |
---|
| 202 | print *,'slab_gm_smax',gm_smax |
---|
| 203 | ! ----------------------------------------------------------- |
---|
| 204 | ! Define ocean / continent mask (no flux into continent cell) |
---|
| 205 | allocate(zmasqu(ip1jmp1)) |
---|
| 206 | allocate(zmasqv(ip1jm)) |
---|
| 207 | zmasqu=1. |
---|
| 208 | zmasqv=1. |
---|
| 209 | |
---|
| 210 | ! convert mask to 2D grid |
---|
| 211 | CALL gr_fi_dyn(1,iip1,jjp1,zmasq,zmasq_2d) |
---|
| 212 | ! put flux mask to 0 at boundaries of continent cells |
---|
| 213 | DO i=1,ip1jmp1-1 |
---|
| 214 | IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+1).GT.1e-5) THEN |
---|
| 215 | zmasqu(i)=0. |
---|
| 216 | ENDIF |
---|
| 217 | END DO |
---|
| 218 | DO i=iip1,ip1jmp1,iip1 |
---|
| 219 | zmasqu(i)=zmasqu(i-iim) |
---|
| 220 | END DO |
---|
| 221 | DO i=1,ip1jm |
---|
| 222 | IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+iip1).GT.1e-5) THEN |
---|
| 223 | zmasqv(i)=0. |
---|
| 224 | END IF |
---|
| 225 | END DO |
---|
| 226 | |
---|
| 227 | ! ----------------------------------------------------------- |
---|
| 228 | ! Coriolis and friction for Ekman transport |
---|
| 229 | ! slab_ekman=2 |
---|
| 230 | slab_ekman=.TRUE. |
---|
| 231 | CALL getin("slab_ekman",slab_ekman) |
---|
| 232 | ! IF (slab_ekman.GT.0) THEN |
---|
| 233 | IF (slab_ekman) THEN |
---|
| 234 | allocate(unsfv(ip1jm)) |
---|
| 235 | allocate(unsev(ip1jm)) |
---|
| 236 | allocate(unsfu(ip1jmp1)) |
---|
| 237 | allocate(unseu(ip1jmp1)) |
---|
| 238 | allocate(unsbv(ip1jm)) |
---|
| 239 | |
---|
| 240 | ! Drage parameter (in s-1) |
---|
| 241 | eps=1e-5 |
---|
| 242 | CALL getin('slab_eps',eps) |
---|
| 243 | print *,'epsilon=',eps |
---|
| 244 | ! Coriolis paraleter |
---|
| 245 | ff=fext*unsairez |
---|
| 246 | |
---|
| 247 | ! coefs to convert tau_x, tau_y to Ekman mass fluxes |
---|
| 248 | ! on 2D grid v points |
---|
| 249 | IF (slab_sverdrup) THEN |
---|
| 250 | ! unsev factor [0 1] for transition Ekman / Sverdrup near the equator (f<<eps) |
---|
| 251 | ! f/eps gives Rossby Radius |
---|
| 252 | DO i=1,ip1jm |
---|
| 253 | unsev(i)=exp(-ff(i)*ff(i)/eps**2) |
---|
| 254 | ENDDO |
---|
| 255 | ELSE |
---|
| 256 | ! Drag vs Coriolis balance |
---|
| 257 | DO i=1,ip1jm |
---|
| 258 | unsev(i)=eps**2/(ff(i)*ff(i)+eps**2) |
---|
| 259 | ENDDO |
---|
| 260 | END IF ! slab_sverdrup |
---|
| 261 | |
---|
| 262 | ! beta(:)=2*omeg*cos(rlatv(:))/rad |
---|
| 263 | !!!!! Here, beta contains rlatv from ini_slab_transp_geom !!!!! |
---|
| 264 | beta(:)=2*omeg*cos(beta(:))/rad |
---|
| 265 | !!!!! |
---|
| 266 | ! 1/beta (derivative of Coriolos parameter) |
---|
| 267 | DO i=1,jjp1-1 |
---|
| 268 | unsbv((i-1)*iip1+1:i*iip1)=unsev((i-1)*iip1+1:i*iip1)/beta(i) |
---|
| 269 | END DO |
---|
| 270 | ! 1/f |
---|
| 271 | ff=SIGN(MAX(ABS(ff),eps/100.),ff) ! avoid value 0 at equator... |
---|
| 272 | DO i=1,ip1jm |
---|
| 273 | unsfv(i)=(1.-unsev(i))/ff(i) |
---|
| 274 | END DO |
---|
| 275 | ! compute values on 2D u grid |
---|
| 276 | ! 1/eps |
---|
| 277 | unsev(:)=unsev(:)/eps |
---|
| 278 | CALL gr_v_scal(1,unsfv,unsfu) |
---|
| 279 | CALL gr_v_scal(1,unsev,unseu) |
---|
| 280 | END IF |
---|
| 281 | |
---|
| 282 | END SUBROUTINE ini_slab_transp |
---|
| 283 | |
---|
| 284 | SUBROUTINE divgrad_phy(nlevs,temp,delta) |
---|
| 285 | ! Computes temperature tendency due to horizontal diffusion : |
---|
| 286 | ! T Laplacian, later multiplied by diffusion coef and time-step |
---|
| 287 | |
---|
| 288 | IMPLICIT NONE |
---|
| 289 | |
---|
| 290 | INTEGER, INTENT(IN) :: nlevs ! nlevs : slab layers |
---|
| 291 | REAL, INTENT(IN) :: temp(klon_glo,nlevs) ! slab temperature |
---|
| 292 | REAL , INTENT(OUT) :: delta(klon_glo,nlevs) ! temp laplacian (heat flux div.) |
---|
| 293 | REAL :: delta_2d((nbp_lon+1)*nbp_lat,nlevs) |
---|
| 294 | REAL ghx((nbp_lon+1)*nbp_lat,nlevs), ghy((nbp_lon+1)*(nbp_lat-1),nlevs) |
---|
| 295 | INTEGER :: ll,iip1,jjp1 |
---|
| 296 | |
---|
| 297 | iip1=nbp_lon+1 |
---|
| 298 | jjp1=nbp_lat |
---|
| 299 | |
---|
| 300 | ! transpose temp to 2D horiz. grid |
---|
| 301 | CALL gr_fi_dyn(nlevs,iip1,jjp1,temp,delta_2d) |
---|
| 302 | ! computes gradient (proportional to heat flx) |
---|
| 303 | CALL grad(nlevs,delta_2d,ghx,ghy) |
---|
| 304 | ! put flux to 0 at ocean / continent boundary |
---|
| 305 | DO ll=1,nlevs |
---|
| 306 | ghx(:,ll)=ghx(:,ll)*zmasqu |
---|
| 307 | ghy(:,ll)=ghy(:,ll)*zmasqv |
---|
| 308 | END DO |
---|
| 309 | ! flux divergence |
---|
| 310 | CALL diverg(nlevs,ghx,ghy,delta_2d) |
---|
| 311 | ! laplacian back to 1D grid |
---|
| 312 | CALL gr_dyn_fi(nlevs,iip1,jjp1,delta_2d,delta) |
---|
| 313 | |
---|
| 314 | RETURN |
---|
| 315 | END SUBROUTINE divgrad_phy |
---|
| 316 | |
---|
| 317 | |
---|
| 318 | SUBROUTINE slab_ekman2(tx_phy,ty_phy,ts_phy,dt_phy_ek,dt_phy_gm,slab_gm) |
---|
| 319 | ! Temperature tendency for 2-layers slab ocean |
---|
| 320 | ! note : tendency dt later multiplied by (delta time)/(rho.H) |
---|
| 321 | ! to convert from divergence of heat fluxes to T |
---|
| 322 | |
---|
| 323 | ! 11/16 : Inclusion of GM-like eddy advection |
---|
| 324 | |
---|
| 325 | IMPLICIT NONE |
---|
| 326 | |
---|
| 327 | LOGICAL,INTENT(in) :: slab_gm |
---|
| 328 | ! Here, temperature and flux variables are on 2 layers |
---|
| 329 | INTEGER ij |
---|
| 330 | |
---|
| 331 | ! wind stress variables |
---|
| 332 | REAL tx_phy(klon_glo),ty_phy(klon_glo) |
---|
| 333 | REAL txv((nbp_lon+1)*(nbp_lat-1)), tyv((nbp_lon+1)*(nbp_lat-1)) |
---|
| 334 | REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat) |
---|
| 335 | REAL tcurl((nbp_lon+1)*(nbp_lat-1)) |
---|
| 336 | ! slab temperature on 1D, 2D grid |
---|
| 337 | REAL ts_phy(klon_glo,2), ts((nbp_lon+1)*nbp_lat,2) |
---|
| 338 | ! Temperature gradient, v-points |
---|
| 339 | REAL dty((nbp_lon+1)*(nbp_lat-1)),dtx((nbp_lon+1)*nbp_lat) |
---|
| 340 | ! Vertical temperature difference, V-points |
---|
| 341 | REAL dtz((nbp_lon+1)*(nbp_lat-1)) |
---|
| 342 | ! zonal and meridional mass fluxes at u, v points (2D grid) |
---|
| 343 | REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1)) |
---|
| 344 | ! vertical mass flux between the 2 layers |
---|
| 345 | REAL fluxv_ek((nbp_lon+1)*nbp_lat) |
---|
| 346 | REAL fluxv_gm((nbp_lon+1)*nbp_lat) |
---|
| 347 | ! zonal and meridional heat fluxes |
---|
| 348 | REAL fluxtz((nbp_lon+1)*nbp_lat,2) |
---|
| 349 | REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2) |
---|
| 350 | ! temperature tendency (in K.s-1.kg.m-2) |
---|
| 351 | REAL dt_ek((nbp_lon+1)*nbp_lat,2), dt_phy_ek(klon_glo,2) |
---|
| 352 | REAL dt_gm((nbp_lon+1)*nbp_lat,2), dt_phy_gm(klon_glo,2) |
---|
| 353 | ! helper vars |
---|
| 354 | REAL zonavg, fluxv |
---|
| 355 | REAL, PARAMETER :: sea_den=1025. ! sea water density |
---|
| 356 | |
---|
| 357 | INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1 |
---|
| 358 | |
---|
| 359 | ! Grid definitions |
---|
| 360 | iim=nbp_lon |
---|
| 361 | iip1=nbp_lon+1 |
---|
| 362 | iip2=nbp_lon+2 |
---|
| 363 | jjp1=nbp_lat |
---|
| 364 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
---|
| 365 | ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1 |
---|
| 366 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
---|
| 367 | ! Convert temperature to 2D grid |
---|
| 368 | CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts) |
---|
| 369 | |
---|
| 370 | ! ------------------------------------ |
---|
| 371 | ! Ekman mass fluxes and Temp tendency |
---|
| 372 | ! ------------------------------------ |
---|
| 373 | ! Convert taux,y to 2D scalar grid |
---|
| 374 | ! north and south poles tx,ty no meaning |
---|
| 375 | tx_phy(1)=0. |
---|
| 376 | tx_phy(klon_glo)=0. |
---|
| 377 | ty_phy(1)=0. |
---|
| 378 | ty_phy(klon_glo)=0. |
---|
| 379 | CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu) |
---|
| 380 | CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu) |
---|
| 381 | |
---|
| 382 | IF (ekman_zonavg) THEN ! use zonal average of wind stress |
---|
| 383 | DO ij=1,jjp1-2 |
---|
| 384 | zonavg=SUM(txu(ij*iip1+1:ij*iip1+iim))/iim |
---|
| 385 | txu(ij*iip1+1:(ij+1)*iip1)=zonavg |
---|
| 386 | zonavg=SUM(tyu(ij*iip1+1:ij*iip1+iim))/iim |
---|
| 387 | tyu(ij*iip1+1:(ij+1)*iip1)=zonavg |
---|
| 388 | END DO |
---|
| 389 | END IF |
---|
| 390 | |
---|
| 391 | ! ---- Meridional flux ---- |
---|
| 392 | ! wind stress at v points |
---|
| 393 | CALL gr_scal_v(1,txu,txv) |
---|
| 394 | |
---|
| 395 | ! Ekman balance |
---|
| 396 | fluxm=-txv*unsfv ! in kg.s-1.m-1 (zonal distance) |
---|
| 397 | |
---|
| 398 | IF (slab_sverdrup) THEN ! Sverdrup bal. near equator |
---|
| 399 | tcurl=(txu(1:ip1jm)-txu(iip2:ip1jmp1))/cv(:) ! dtx/dy |
---|
| 400 | !poles curl = 0 |
---|
| 401 | tcurl(1:iip1)=0. |
---|
| 402 | tcurl(ip1jmi1+1:ip1jm)=0. |
---|
| 403 | fluxm=fluxm-tcurl*unsbv |
---|
| 404 | ENDIF |
---|
| 405 | |
---|
| 406 | IF (slab_tyeq) THEN ! frictional meridional transport at equator |
---|
| 407 | CALL gr_scal_v(1,tyu,tyv) |
---|
| 408 | fluxm=fluxm+tyv*unsev ! in kg.s-1.m-1 (zonal distance) |
---|
| 409 | ENDIF |
---|
| 410 | |
---|
| 411 | ! apply continent mask, multiply by horiz grid dimension |
---|
| 412 | fluxm=fluxm*cv*cuvsurcv*zmasqv |
---|
| 413 | |
---|
| 414 | ! ---- Zonal flux ---- |
---|
| 415 | IF (ekman_zonadv) THEN |
---|
| 416 | CALL gr_scal_u(1,txu,txu) ! wind stress at u points |
---|
| 417 | CALL gr_scal_u(1,tyu,tyu) |
---|
| 418 | fluxz=tyu*unsfu+txu*unseu ! Ekman + frictional flow |
---|
| 419 | ! apply continent mask, multiply by horiz grid dimension |
---|
| 420 | fluxz=fluxz*cu*cvusurcu*zmasqu |
---|
| 421 | END IF |
---|
| 422 | |
---|
| 423 | ! ---- Vertical mass flux ---- |
---|
| 424 | ! from mass budget (divergence of horiz fluxes) |
---|
| 425 | IF (ekman_zonadv) THEN |
---|
| 426 | DO ij=iip2,ip1jm |
---|
| 427 | fluxv_ek(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1) |
---|
| 428 | ENDDO |
---|
| 429 | ELSE |
---|
| 430 | DO ij=iip2,ip1jm |
---|
| 431 | fluxv_ek(ij)=-fluxm(ij)+fluxm(ij-iip1) |
---|
| 432 | ENDDO |
---|
| 433 | END IF |
---|
| 434 | DO ij=iip1,ip1jmi1,iip1 |
---|
| 435 | fluxv_ek(ij+1)=fluxv_ek(ij+iip1) |
---|
| 436 | END DO |
---|
| 437 | ! vertical mass flux at Poles |
---|
| 438 | fluxv_ek(1)=-SUM(fluxm(1:iim)) |
---|
| 439 | fluxv_ek(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1)) |
---|
| 440 | |
---|
| 441 | ! ---- Meridional heat fluxes ---- |
---|
| 442 | DO ij=1,ip1jm |
---|
| 443 | ! centered scheme |
---|
| 444 | fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2. |
---|
| 445 | fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2. |
---|
| 446 | END DO |
---|
| 447 | |
---|
| 448 | ! ---- Zonal heat fluxes ---- |
---|
| 449 | ! Schema upstream |
---|
| 450 | IF (ekman_zonadv) THEN |
---|
| 451 | fluxtz(:,:)=0. |
---|
| 452 | DO ij=iip2,ip1jm |
---|
| 453 | IF (fluxz(ij).GE.0.) THEN |
---|
| 454 | fluxtz(ij,1)=fluxz(ij)*ts(ij,1) |
---|
| 455 | fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2) |
---|
| 456 | ELSE |
---|
| 457 | fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1) |
---|
| 458 | fluxtz(ij,2)=-fluxz(ij)*ts(ij,2) |
---|
| 459 | ENDIF |
---|
| 460 | ENDDO |
---|
| 461 | DO ij=iip1,ip1jmp1,iip1 |
---|
| 462 | fluxtz(ij,:)=fluxtz(ij-iim,:) |
---|
| 463 | END DO |
---|
| 464 | ELSE |
---|
| 465 | fluxtz(:,:)=0. |
---|
| 466 | ENDIF |
---|
| 467 | |
---|
| 468 | ! ---- Temperature tendency, horizontal advection: ---- |
---|
| 469 | DO ij=iip2,ip1jm |
---|
| 470 | dt_ek(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) & |
---|
| 471 | +fluxtm(ij,:)-fluxtm(ij-iip1,:) |
---|
| 472 | END DO |
---|
| 473 | ! Poles |
---|
| 474 | dt_ek(1,:)=SUM(fluxtm(1:iim,:),dim=1) |
---|
| 475 | dt_ek(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1) |
---|
| 476 | |
---|
| 477 | ! ------------------------------------ |
---|
| 478 | ! GM mass fluxes and Temp tendency |
---|
| 479 | ! ------------------------------------ |
---|
| 480 | IF (slab_gm) THEN |
---|
| 481 | ! Vertical Temperature difference T1-T2 on v-grid points |
---|
| 482 | CALL gr_scal_v(1,ts(:,1)-ts(:,2),dtz) |
---|
| 483 | dtz(:)=MAX(dtz(:),0.25) |
---|
| 484 | ! Horizontal Temperature differences |
---|
| 485 | CALL grad(1,(ts(:,1)+ts(:,2))/2.,dtx,dty) |
---|
| 486 | ! Meridional flux = -k.s (s=dyT/dzT) |
---|
| 487 | ! Continent mask, multiply by dz/dy |
---|
| 488 | fluxm=dty/dtz*500.*cuvsurcv*zmasqv |
---|
| 489 | ! slope limitation, multiply by kappa |
---|
| 490 | fluxm=-gmkappa*SIGN(MIN(ABS(fluxm),gm_smax*cv*cuvsurcv),dty) |
---|
| 491 | ! convert to kg/s |
---|
| 492 | fluxm(:)=fluxm(:)*sea_den |
---|
| 493 | ! Zonal flux = 0. (temporary) |
---|
| 494 | fluxz(:)=0. |
---|
| 495 | ! Vertical mass flux from mass budget (divergence of horiz fluxes) |
---|
| 496 | DO ij=iip2,ip1jm |
---|
| 497 | fluxv_gm(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1) |
---|
| 498 | ENDDO |
---|
| 499 | DO ij=iip1,ip1jmi1,iip1 |
---|
| 500 | fluxv_gm(ij+1)=fluxv_gm(ij+iip1) |
---|
| 501 | END DO |
---|
| 502 | ! vertical mass flux at Poles |
---|
| 503 | fluxv_gm(1)=-SUM(fluxm(1:iim)) |
---|
| 504 | fluxv_gm(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1)) |
---|
| 505 | |
---|
| 506 | ! Meridional heat fluxes |
---|
| 507 | DO ij=1,ip1jm |
---|
| 508 | ! centered scheme |
---|
| 509 | fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2. |
---|
| 510 | fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2. |
---|
| 511 | END DO |
---|
| 512 | |
---|
| 513 | ! Zonal heat fluxes |
---|
| 514 | ! Schema upstream |
---|
| 515 | DO ij=iip2,ip1jm |
---|
| 516 | IF (fluxz(ij).GE.0.) THEN |
---|
| 517 | fluxtz(ij,1)=fluxz(ij)*ts(ij,1) |
---|
| 518 | fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2) |
---|
| 519 | ELSE |
---|
| 520 | fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1) |
---|
| 521 | fluxtz(ij,2)=-fluxz(ij)*ts(ij,2) |
---|
| 522 | ENDIF |
---|
| 523 | ENDDO |
---|
| 524 | DO ij=iip1*2,ip1jmp1,iip1 |
---|
| 525 | fluxtz(ij,:)=fluxtz(ij-iim,:) |
---|
| 526 | END DO |
---|
| 527 | |
---|
| 528 | ! Temperature tendency : |
---|
| 529 | ! divergence of horizontal heat fluxes |
---|
| 530 | DO ij=iip2,ip1jm |
---|
| 531 | dt_gm(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) & |
---|
| 532 | +fluxtm(ij,:)-fluxtm(ij-iip1,:) |
---|
| 533 | END DO |
---|
| 534 | ! Poles |
---|
| 535 | dt_gm(1,:)=SUM(fluxtm(1:iim,:),dim=1) |
---|
| 536 | dt_gm(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1) |
---|
| 537 | ELSE |
---|
| 538 | dt_gm(:,:)=0. |
---|
| 539 | fluxv_gm(:)=0. |
---|
| 540 | ENDIF ! slab_gm |
---|
| 541 | |
---|
| 542 | ! ------------------------------------ |
---|
| 543 | ! Temp tendency from vertical advection |
---|
| 544 | ! Divide by cell area |
---|
| 545 | ! ------------------------------------ |
---|
| 546 | ! vertical heat flux = mass flux * T, upstream scheme |
---|
| 547 | DO ij=iip2,ip1jm |
---|
| 548 | ! net flux, needed for upstream scheme |
---|
| 549 | fluxv=fluxv_ek(ij)+fluxv_gm(ij) |
---|
| 550 | ! temperature tendencies from vertical fluxes |
---|
| 551 | IF (fluxv.GT.0.) THEN |
---|
| 552 | dt_ek(ij,1)=dt_ek(ij,1)+fluxv_ek(ij)*ts(ij,2) |
---|
| 553 | dt_ek(ij,2)=dt_ek(ij,2)-fluxv_ek(ij)*ts(ij,2) |
---|
| 554 | dt_gm(ij,1)=dt_gm(ij,1)+fluxv_gm(ij)*ts(ij,2) |
---|
| 555 | dt_gm(ij,2)=dt_gm(ij,2)-fluxv_gm(ij)*ts(ij,2) |
---|
| 556 | ELSE |
---|
| 557 | dt_ek(ij,1)=dt_ek(ij,1)+fluxv_ek(ij)*ts(ij,1) |
---|
| 558 | dt_ek(ij,2)=dt_ek(ij,2)-fluxv_ek(ij)*ts(ij,1) |
---|
| 559 | dt_gm(ij,1)=dt_gm(ij,1)+fluxv_gm(ij)*ts(ij,1) |
---|
| 560 | dt_gm(ij,2)=dt_gm(ij,2)-fluxv_gm(ij)*ts(ij,1) |
---|
| 561 | ENDIF |
---|
| 562 | ! divide by cell area |
---|
| 563 | dt_ek(ij,:)=dt_ek(ij,:)/aire(ij) |
---|
| 564 | dt_gm(ij,:)=dt_gm(ij,:)/aire(ij) |
---|
| 565 | END DO |
---|
| 566 | ! North Pole |
---|
| 567 | fluxv=fluxv_ek(1)+fluxv_gm(1) |
---|
| 568 | IF (fluxv.GT.0.) THEN |
---|
| 569 | dt_ek(1,1)=dt_ek(1,1)+fluxv_ek(1)*ts(1,2) |
---|
| 570 | dt_ek(1,2)=dt_ek(1,2)-fluxv_ek(1)*ts(1,2) |
---|
| 571 | dt_gm(1,1)=dt_gm(1,1)+fluxv_gm(1)*ts(1,2) |
---|
| 572 | dt_gm(1,2)=dt_gm(1,2)-fluxv_gm(1)*ts(1,2) |
---|
| 573 | ELSE |
---|
| 574 | dt_ek(1,1)=dt_ek(1,1)+fluxv_ek(1)*ts(1,1) |
---|
| 575 | dt_ek(1,2)=dt_ek(1,2)-fluxv_ek(1)*ts(1,1) |
---|
| 576 | dt_gm(1,1)=dt_gm(1,1)+fluxv_gm(1)*ts(1,1) |
---|
| 577 | dt_gm(1,2)=dt_gm(1,2)-fluxv_gm(1)*ts(1,1) |
---|
| 578 | ENDIF |
---|
| 579 | dt_ek(1,:)=dt_ek(1,:)/apoln |
---|
| 580 | dt_gm(1,:)=dt_gm(1,:)/apoln |
---|
| 581 | ! South pole |
---|
| 582 | fluxv=fluxv_ek(ip1jmp1)+fluxv_gm(ip1jmp1) |
---|
| 583 | IF (fluxv.GT.0.) THEN |
---|
| 584 | dt_ek(ip1jmp1,1)=dt_ek(ip1jmp1,1)+fluxv_ek(ip1jmp1)*ts(ip1jmp1,2) |
---|
| 585 | dt_ek(ip1jmp1,2)=dt_ek(ip1jmp1,2)-fluxv_ek(ip1jmp1)*ts(ip1jmp1,2) |
---|
| 586 | dt_gm(ip1jmp1,1)=dt_gm(ip1jmp1,1)+fluxv_gm(ip1jmp1)*ts(ip1jmp1,2) |
---|
| 587 | dt_gm(ip1jmp1,2)=dt_gm(ip1jmp1,2)-fluxv_gm(ip1jmp1)*ts(ip1jmp1,2) |
---|
| 588 | ELSE |
---|
| 589 | dt_ek(ip1jmp1,1)=dt_ek(ip1jmp1,1)+fluxv_ek(ip1jmp1)*ts(ip1jmp1,1) |
---|
| 590 | dt_ek(ip1jmp1,2)=dt_ek(ip1jmp1,2)-fluxv_ek(ip1jmp1)*ts(ip1jmp1,1) |
---|
| 591 | dt_gm(ip1jmp1,1)=dt_gm(ip1jmp1,1)+fluxv_gm(ip1jmp1)*ts(ip1jmp1,1) |
---|
| 592 | dt_gm(ip1jmp1,2)=dt_gm(ip1jmp1,2)-fluxv_gm(ip1jmp1)*ts(ip1jmp1,1) |
---|
| 593 | ENDIF |
---|
| 594 | dt_ek(ip1jmp1,:)=dt_ek(ip1jmp1,:)/apols |
---|
| 595 | dt_gm(ip1jmp1,:)=dt_gm(ip1jmp1,:)/apols |
---|
| 596 | |
---|
| 597 | dt_ek(2:iip1,1)=dt_ek(1,1) |
---|
| 598 | dt_ek(2:iip1,2)=dt_ek(1,2) |
---|
| 599 | dt_gm(2:iip1,1)=dt_gm(1,1) |
---|
| 600 | dt_gm(2:iip1,2)=dt_gm(1,2) |
---|
| 601 | dt_ek(ip1jm+1:ip1jmp1,1)=dt_ek(ip1jmp1,1) |
---|
| 602 | dt_ek(ip1jm+1:ip1jmp1,2)=dt_ek(ip1jmp1,2) |
---|
| 603 | dt_gm(ip1jm+1:ip1jmp1,1)=dt_gm(ip1jmp1,1) |
---|
| 604 | dt_gm(ip1jm+1:ip1jmp1,2)=dt_gm(ip1jmp1,2) |
---|
| 605 | |
---|
| 606 | DO ij=iip1,ip1jmi1,iip1 |
---|
| 607 | dt_gm(ij+1,:)=dt_gm(ij+iip1,:) |
---|
| 608 | dt_ek(ij+1,:)=dt_ek(ij+iip1,:) |
---|
| 609 | END DO |
---|
| 610 | |
---|
| 611 | ! T tendency back to 1D grid... |
---|
| 612 | CALL gr_dyn_fi(2,iip1,jjp1,dt_ek,dt_phy_ek) |
---|
| 613 | CALL gr_dyn_fi(2,iip1,jjp1,dt_gm,dt_phy_gm) |
---|
| 614 | |
---|
| 615 | RETURN |
---|
| 616 | END SUBROUTINE slab_ekman2 |
---|
| 617 | |
---|
| 618 | SUBROUTINE slab_gmdiff(ts_phy,dt_phy) |
---|
| 619 | ! Temperature tendency for 2-layers slab ocean |
---|
| 620 | ! Due to Gent-McWilliams type eddy-induced advection |
---|
| 621 | ! Not used anymore ? (included in Ekman) |
---|
| 622 | |
---|
| 623 | IMPLICIT NONE |
---|
| 624 | |
---|
| 625 | ! Here, temperature and flux variables are on 2 layers |
---|
| 626 | INTEGER ij |
---|
| 627 | ! Temperature gradient, v-points |
---|
| 628 | REAL dty((nbp_lon+1)*(nbp_lat-1)),dtx((nbp_lon+1)*nbp_lat) |
---|
| 629 | ! Vertical temperature difference, V-points |
---|
| 630 | REAL dtz((nbp_lon+1)*(nbp_lat-1)) |
---|
| 631 | ! slab temperature on 1D, 2D grid |
---|
| 632 | REAL ts_phy(klon_glo,2),ts((nbp_lon+1)*nbp_lat,2) |
---|
| 633 | ! zonal and meridional mass fluxes at u, v points (2D grid) |
---|
| 634 | REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1)) |
---|
| 635 | ! vertical mass flux between the 2 layers |
---|
| 636 | REAL fluxv((nbp_lon+1)*nbp_lat) |
---|
| 637 | ! zonal and meridional heat fluxes |
---|
| 638 | REAL fluxtz((nbp_lon+1)*nbp_lat,2) |
---|
| 639 | REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2) |
---|
| 640 | ! temperature tendency (in K.s-1.kg.m-2) |
---|
| 641 | REAL dt((nbp_lon+1)*nbp_lat,2), dt_phy(klon_glo,2) |
---|
| 642 | |
---|
| 643 | INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1 |
---|
| 644 | |
---|
| 645 | ! Grid definitions |
---|
| 646 | iim=nbp_lon |
---|
| 647 | iip1=nbp_lon+1 |
---|
| 648 | iip2=nbp_lon+2 |
---|
| 649 | jjp1=nbp_lat |
---|
| 650 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
---|
| 651 | ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1 |
---|
| 652 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
---|
| 653 | |
---|
| 654 | ! Convert temperature to 2D grid |
---|
| 655 | CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts) |
---|
| 656 | ! Vertical Temperature difference T1-T2 on v-grid points |
---|
| 657 | CALL gr_scal_v(1,ts(:,1)-ts(:,2),dtz) |
---|
| 658 | dtz(:)=MAX(dtz(:),0.25) |
---|
| 659 | ! Horizontal Temperature differences |
---|
| 660 | CALL grad(1,(ts(:,1)+ts(:,2))/2.,dtx,dty) |
---|
| 661 | ! Meridional flux = -k.s (s=dyT/dzT) |
---|
| 662 | ! Continent mask, multiply by dz/dy |
---|
| 663 | fluxm=dty/dtz*500.*cuvsurcv*zmasqv |
---|
| 664 | ! slope limitation, multiply by kappa |
---|
| 665 | fluxm=-gmkappa*SIGN(MIN(ABS(fluxm),gm_smax*cv*cuvsurcv),dty) |
---|
| 666 | ! Zonal flux = 0. (temporary) |
---|
| 667 | fluxz(:)=0. |
---|
| 668 | ! Vertical mass flux from mass budget (divergence of horiz fluxes) |
---|
| 669 | DO ij=iip2,ip1jm |
---|
| 670 | fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1) |
---|
| 671 | ENDDO |
---|
| 672 | DO ij=iip1,ip1jmi1,iip1 |
---|
| 673 | fluxv(ij+1)=fluxv(ij+iip1) |
---|
| 674 | END DO |
---|
| 675 | ! vertical mass flux at Poles |
---|
| 676 | fluxv(1)=-SUM(fluxm(1:iim)) |
---|
| 677 | fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1)) |
---|
| 678 | fluxv=fluxv |
---|
| 679 | |
---|
| 680 | ! Meridional heat fluxes |
---|
| 681 | DO ij=1,ip1jm |
---|
| 682 | ! centered scheme |
---|
| 683 | fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2. |
---|
| 684 | fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2. |
---|
| 685 | END DO |
---|
| 686 | |
---|
| 687 | ! Zonal heat fluxes |
---|
| 688 | ! Schema upstream |
---|
| 689 | DO ij=iip2,ip1jm |
---|
| 690 | IF (fluxz(ij).GE.0.) THEN |
---|
| 691 | fluxtz(ij,1)=fluxz(ij)*ts(ij,1) |
---|
| 692 | fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2) |
---|
| 693 | ELSE |
---|
| 694 | fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1) |
---|
| 695 | fluxtz(ij,2)=-fluxz(ij)*ts(ij,2) |
---|
| 696 | ENDIF |
---|
| 697 | ENDDO |
---|
| 698 | DO ij=iip1*2,ip1jmp1,iip1 |
---|
| 699 | fluxtz(ij,:)=fluxtz(ij-iim,:) |
---|
| 700 | END DO |
---|
| 701 | |
---|
| 702 | ! Temperature tendency : |
---|
| 703 | DO ij=iip2,ip1jm |
---|
| 704 | ! divergence of horizontal heat fluxes |
---|
| 705 | dt(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) & |
---|
| 706 | +fluxtm(ij,:)-fluxtm(ij-iip1,:) |
---|
| 707 | ! + vertical heat flux (mass flux * T, upstream scheme) |
---|
| 708 | IF (fluxv(ij).GT.0.) THEN |
---|
| 709 | dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,2) |
---|
| 710 | dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,2) |
---|
| 711 | ELSE |
---|
| 712 | dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,1) |
---|
| 713 | dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,1) |
---|
| 714 | ENDIF |
---|
| 715 | ! divide by cell area |
---|
| 716 | dt(ij,:)=dt(ij,:)/aire(ij) |
---|
| 717 | END DO |
---|
| 718 | DO ij=iip1,ip1jmi1,iip1 |
---|
| 719 | dt(ij+1,:)=dt(ij+iip1,:) |
---|
| 720 | END DO |
---|
| 721 | ! Poles |
---|
| 722 | dt(1,:)=SUM(fluxtm(1:iim,:),dim=1) |
---|
| 723 | IF (fluxv(1).GT.0.) THEN |
---|
| 724 | dt(1,1)=dt(1,1)+fluxv(1)*ts(1,2) |
---|
| 725 | dt(1,2)=dt(1,2)-fluxv(1)*ts(1,2) |
---|
| 726 | ELSE |
---|
| 727 | dt(1,1)=dt(1,1)+fluxv(1)*ts(1,1) |
---|
| 728 | dt(1,2)=dt(1,2)-fluxv(1)*ts(1,1) |
---|
| 729 | ENDIF |
---|
| 730 | dt(1,:)=dt(1,:)/apoln |
---|
| 731 | dt(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1) |
---|
| 732 | IF (fluxv(ip1jmp1).GT.0.) THEN |
---|
| 733 | dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,2) |
---|
| 734 | dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,2) |
---|
| 735 | ELSE |
---|
| 736 | dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,1) |
---|
| 737 | dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,1) |
---|
| 738 | ENDIF |
---|
| 739 | dt(ip1jmp1,:)=dt(ip1jmp1,:)/apols |
---|
| 740 | dt(2:iip1,1)=dt(1,1) |
---|
| 741 | dt(2:iip1,2)=dt(1,2) |
---|
| 742 | dt(ip1jm+1:ip1jmp1,1)=dt(ip1jmp1,1) |
---|
| 743 | dt(ip1jm+1:ip1jmp1,2)=dt(ip1jmp1,2) |
---|
| 744 | |
---|
| 745 | ! T tendency back to 1D grid... |
---|
| 746 | CALL gr_dyn_fi(2,iip1,jjp1,dt,dt_phy) |
---|
| 747 | |
---|
| 748 | RETURN |
---|
| 749 | END SUBROUTINE slab_gmdiff |
---|
| 750 | |
---|
| 751 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 752 | !!! Routines copied from dynamics !!! |
---|
| 753 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 754 | |
---|
| 755 | SUBROUTINE gr_fi_dyn(nfield,im,jm,pfi,pdyn) |
---|
| 756 | ! Transfer a variable from 1D "physics" grid to 2D "dynamics" grid |
---|
| 757 | IMPLICIT NONE |
---|
| 758 | |
---|
| 759 | INTEGER,INTENT(IN) :: im,jm,nfield |
---|
| 760 | REAL,INTENT(IN) :: pfi(klon_glo,nfield) ! on 1D grid |
---|
| 761 | REAL,INTENT(OUT) :: pdyn(im,jm,nfield) ! on 2D grid |
---|
| 762 | |
---|
| 763 | INTEGER :: i,j,ifield,ig |
---|
| 764 | |
---|
| 765 | DO ifield=1,nfield |
---|
| 766 | ! Handle poles |
---|
| 767 | DO i=1,im |
---|
| 768 | pdyn(i,1,ifield)=pfi(1,ifield) |
---|
| 769 | pdyn(i,jm,ifield)=pfi(klon_glo,ifield) |
---|
| 770 | ENDDO |
---|
| 771 | ! Other points |
---|
| 772 | DO j=2,jm-1 |
---|
| 773 | ig=2+(j-2)*(im-1) |
---|
| 774 | CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1) |
---|
| 775 | pdyn(im,j,ifield)=pdyn(1,j,ifield) |
---|
| 776 | ENDDO |
---|
| 777 | ENDDO ! of DO ifield=1,nfield |
---|
| 778 | |
---|
| 779 | END SUBROUTINE gr_fi_dyn |
---|
| 780 | |
---|
| 781 | SUBROUTINE gr_dyn_fi(nfield,im,jm,pdyn,pfi) |
---|
| 782 | ! Transfer a variable from 2D "dynamics" grid to 1D "physics" grid |
---|
| 783 | IMPLICIT NONE |
---|
| 784 | |
---|
| 785 | INTEGER,INTENT(IN) :: im,jm,nfield |
---|
| 786 | REAL,INTENT(IN) :: pdyn(im,jm,nfield) ! on 2D grid |
---|
| 787 | REAL,INTENT(OUT) :: pfi(klon_glo,nfield) ! on 1D grid |
---|
| 788 | |
---|
| 789 | INTEGER j,ifield,ig |
---|
| 790 | |
---|
| 791 | ! Sanity check: |
---|
| 792 | IF(klon_glo.NE.2+(jm-2)*(im-1)) THEN |
---|
| 793 | WRITE(*,*) "gr_dyn_fi error, wrong sizes" |
---|
| 794 | STOP |
---|
| 795 | ENDIF |
---|
| 796 | |
---|
| 797 | ! Handle poles |
---|
| 798 | CALL SCOPY(nfield,pdyn,im*jm,pfi,klon_glo) |
---|
| 799 | CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(klon_glo,1),klon_glo) |
---|
| 800 | ! Other points |
---|
| 801 | DO ifield=1,nfield |
---|
| 802 | DO j=2,jm-1 |
---|
| 803 | ig=2+(j-2)*(im-1) |
---|
| 804 | CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1) |
---|
| 805 | ENDDO |
---|
| 806 | ENDDO |
---|
| 807 | |
---|
| 808 | END SUBROUTINE gr_dyn_fi |
---|
| 809 | |
---|
| 810 | SUBROUTINE grad(klevel,pg,pgx,pgy) |
---|
| 811 | ! compute the covariant components pgx,pgy of the gradient of pg |
---|
| 812 | ! pgx = d(pg)/dx * delta(x) = delta(pg) |
---|
| 813 | IMPLICIT NONE |
---|
| 814 | |
---|
| 815 | INTEGER,INTENT(IN) :: klevel |
---|
| 816 | REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel) |
---|
| 817 | REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel) |
---|
| 818 | REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel) |
---|
| 819 | |
---|
| 820 | INTEGER :: l,ij |
---|
| 821 | INTEGER :: iim,iip1,ip1jm,ip1jmp1 |
---|
| 822 | |
---|
| 823 | iim=nbp_lon |
---|
| 824 | iip1=nbp_lon+1 |
---|
| 825 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
---|
| 826 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
---|
| 827 | |
---|
| 828 | DO l=1,klevel |
---|
| 829 | DO ij=1,ip1jmp1-1 |
---|
| 830 | pgx(ij,l)=pg(ij+1,l)-pg(ij,l) |
---|
| 831 | ENDDO |
---|
| 832 | ! correction for pgx(ip1,j,l) ... |
---|
| 833 | ! ... pgx(iip1,j,l)=pgx(1,j,l) ... |
---|
| 834 | DO ij=iip1,ip1jmp1,iip1 |
---|
| 835 | pgx(ij,l)=pgx(ij-iim,l) |
---|
| 836 | ENDDO |
---|
| 837 | DO ij=1,ip1jm |
---|
| 838 | pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l) |
---|
| 839 | ENDDO |
---|
| 840 | ENDDO |
---|
| 841 | |
---|
| 842 | END SUBROUTINE grad |
---|
| 843 | |
---|
| 844 | SUBROUTINE diverg(klevel,x,y,div) |
---|
| 845 | ! computes the divergence of a vector field of components |
---|
| 846 | ! x,y. x and y being covariant components |
---|
| 847 | IMPLICIT NONE |
---|
| 848 | |
---|
| 849 | INTEGER,INTENT(IN) :: klevel |
---|
| 850 | REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel) |
---|
| 851 | REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel) |
---|
| 852 | REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel) |
---|
| 853 | |
---|
| 854 | INTEGER :: l,ij |
---|
| 855 | INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1 |
---|
| 856 | |
---|
| 857 | REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1) |
---|
| 858 | REAL :: sumypn,sumyps |
---|
| 859 | REAL,EXTERNAL :: SSUM |
---|
| 860 | |
---|
| 861 | iim=nbp_lon |
---|
| 862 | iip1=nbp_lon+1 |
---|
| 863 | iip2=nbp_lon+2 |
---|
| 864 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
---|
| 865 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
---|
| 866 | ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1 |
---|
| 867 | |
---|
| 868 | DO l=1,klevel |
---|
| 869 | DO ij=iip2,ip1jm-1 |
---|
| 870 | div(ij+1,l)= & |
---|
| 871 | cvusurcu(ij+1)*x(ij+1,l)-cvusurcu(ij)*x(ij,l)+ & |
---|
| 872 | cuvsurcv(ij-iim)*y(ij-iim,l)-cuvsurcv(ij+1)*y(ij+1,l) |
---|
| 873 | ENDDO |
---|
| 874 | ! correction for div(1,j,l) ... |
---|
| 875 | ! ... div(1,j,l)= div(iip1,j,l) ... |
---|
| 876 | DO ij=iip2,ip1jm,iip1 |
---|
| 877 | div(ij,l)=div(ij+iim,l) |
---|
| 878 | ENDDO |
---|
| 879 | ! at the poles |
---|
| 880 | DO ij=1,iim |
---|
| 881 | aiy1(ij)=cuvsurcv(ij)*y(ij,l) |
---|
| 882 | aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l) |
---|
| 883 | ENDDO |
---|
| 884 | sumypn=SSUM(iim,aiy1,1)/apoln |
---|
| 885 | sumyps=SSUM(iim,aiy2,1)/apols |
---|
| 886 | DO ij=1,iip1 |
---|
| 887 | div(ij,l)=-sumypn |
---|
| 888 | div(ij+ip1jm,l)=sumyps |
---|
| 889 | ENDDO |
---|
| 890 | ! End (poles) |
---|
| 891 | ENDDO ! of DO l=1,klevel |
---|
| 892 | |
---|
| 893 | !!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 ) |
---|
| 894 | DO l=1,klevel |
---|
| 895 | DO ij=iip2,ip1jm |
---|
| 896 | div(ij,l)=div(ij,l)*unsaire(ij) |
---|
| 897 | ENDDO |
---|
| 898 | ENDDO |
---|
| 899 | |
---|
| 900 | END SUBROUTINE diverg |
---|
| 901 | |
---|
| 902 | SUBROUTINE gr_v_scal(nx,x_v,x_scal) |
---|
| 903 | ! convert values from v points to scalar points on C-grid |
---|
| 904 | ! used to compute unsfu, unseu (u points, but depends only on latitude) |
---|
| 905 | IMPLICIT NONE |
---|
| 906 | |
---|
| 907 | INTEGER,INTENT(IN) :: nx ! number of levels or fields |
---|
| 908 | REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx) |
---|
| 909 | REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx) |
---|
| 910 | |
---|
| 911 | INTEGER :: l,ij |
---|
| 912 | INTEGER :: iip1,iip2,ip1jm,ip1jmp1 |
---|
| 913 | |
---|
| 914 | iip1=nbp_lon+1 |
---|
| 915 | iip2=nbp_lon+2 |
---|
| 916 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
---|
| 917 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
---|
| 918 | |
---|
| 919 | DO l=1,nx |
---|
| 920 | DO ij=iip2,ip1jm |
---|
| 921 | x_scal(ij,l)= & |
---|
| 922 | (airev(ij-iip1)*x_v(ij-iip1,l)+airev(ij)*x_v(ij,l)) & |
---|
| 923 | /(airev(ij-iip1)+airev(ij)) |
---|
| 924 | ENDDO |
---|
| 925 | DO ij=1,iip1 |
---|
| 926 | x_scal(ij,l)=0. |
---|
| 927 | ENDDO |
---|
| 928 | DO ij=ip1jm+1,ip1jmp1 |
---|
| 929 | x_scal(ij,l)=0. |
---|
| 930 | ENDDO |
---|
| 931 | ENDDO |
---|
| 932 | |
---|
| 933 | END SUBROUTINE gr_v_scal |
---|
| 934 | |
---|
| 935 | SUBROUTINE gr_scal_v(nx,x_scal,x_v) |
---|
| 936 | ! convert values from scalar points to v points on C-grid |
---|
| 937 | ! used to compute wind stress at V points |
---|
| 938 | IMPLICIT NONE |
---|
| 939 | |
---|
| 940 | INTEGER,INTENT(IN) :: nx ! number of levels or fields |
---|
| 941 | REAL,INTENT(OUT) :: x_v((nbp_lon+1)*(nbp_lat-1),nx) |
---|
| 942 | REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx) |
---|
| 943 | |
---|
| 944 | INTEGER :: l,ij |
---|
| 945 | INTEGER :: iip1,ip1jm |
---|
| 946 | |
---|
| 947 | iip1=nbp_lon+1 |
---|
| 948 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
---|
| 949 | |
---|
| 950 | DO l=1,nx |
---|
| 951 | DO ij=1,ip1jm |
---|
| 952 | x_v(ij,l)= & |
---|
| 953 | (cu(ij)*cvusurcu(ij)*x_scal(ij,l)+ & |
---|
| 954 | cu(ij+iip1)*cvusurcu(ij+iip1)*x_scal(ij+iip1,l)) & |
---|
| 955 | /(cu(ij)*cvusurcu(ij)+cu(ij+iip1)*cvusurcu(ij+iip1)) |
---|
| 956 | ENDDO |
---|
| 957 | ENDDO |
---|
| 958 | |
---|
| 959 | END SUBROUTINE gr_scal_v |
---|
| 960 | |
---|
| 961 | SUBROUTINE gr_scal_u(nx,x_scal,x_u) |
---|
| 962 | ! convert values from scalar points to U points on C-grid |
---|
| 963 | ! used to compute wind stress at U points |
---|
| 964 | IMPLICIT NONE |
---|
| 965 | |
---|
| 966 | INTEGER,INTENT(IN) :: nx |
---|
| 967 | REAL,INTENT(OUT) :: x_u((nbp_lon+1)*nbp_lat,nx) |
---|
| 968 | REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx) |
---|
| 969 | |
---|
| 970 | INTEGER :: l,ij |
---|
| 971 | INTEGER :: iip1,jjp1,ip1jmp1 |
---|
| 972 | |
---|
| 973 | iip1=nbp_lon+1 |
---|
| 974 | jjp1=nbp_lat |
---|
| 975 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
---|
| 976 | |
---|
| 977 | DO l=1,nx |
---|
| 978 | DO ij=1,ip1jmp1-1 |
---|
| 979 | x_u(ij,l)= & |
---|
| 980 | (aire(ij)*x_scal(ij,l)+aire(ij+1)*x_scal(ij+1,l)) & |
---|
| 981 | /(aire(ij)+aire(ij+1)) |
---|
| 982 | ENDDO |
---|
| 983 | ENDDO |
---|
| 984 | |
---|
| 985 | CALL SCOPY(nx*jjp1,x_u(1,1),iip1,x_u(iip1,1),iip1) |
---|
| 986 | |
---|
| 987 | END SUBROUTINE gr_scal_u |
---|
| 988 | |
---|
| 989 | END MODULE slab_heat_transp_mod |
---|