| 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 :: unsairez(:) ! 1/cell area |
|---|
| 14 | !$OMP THREADPRIVATE(unsairez) |
|---|
| 15 | REAL,SAVE,ALLOCATABLE :: unsaire(:) |
|---|
| 16 | !$OMP THREADPRIVATE(unsaire) |
|---|
| 17 | REAL,SAVE,ALLOCATABLE :: cu(:) ! cell longitude dim (m) |
|---|
| 18 | !$OMP THREADPRIVATE(cu) |
|---|
| 19 | REAL,SAVE,ALLOCATABLE :: cv(:) ! cell latitude dim (m) |
|---|
| 20 | !$OMP THREADPRIVATE(cv) |
|---|
| 21 | REAL,SAVE,ALLOCATABLE :: cuvsurcv(:) ! cu/cv (v points) |
|---|
| 22 | !$OMP THREADPRIVATE(cuvsurcv) |
|---|
| 23 | REAL,SAVE,ALLOCATABLE :: cvusurcu(:) ! cv/cu (u points) |
|---|
| 24 | !$OMP THREADPRIVATE(cvusurcu) |
|---|
| 25 | REAL,SAVE,ALLOCATABLE :: aire(:) ! cell area |
|---|
| 26 | !$OMP THREADPRIVATE(aire) |
|---|
| 27 | REAL,SAVE :: apoln ! area of north pole points |
|---|
| 28 | !$OMP THREADPRIVATE(apoln) |
|---|
| 29 | REAL,SAVE :: apols ! area of south pole points |
|---|
| 30 | !$OMP THREADPRIVATE(apols) |
|---|
| 31 | REAL,SAVE,ALLOCATABLE :: aireu(:) ! area of u cells |
|---|
| 32 | !$OMP THREADPRIVATE(aireu) |
|---|
| 33 | REAL,SAVE,ALLOCATABLE :: airev(:) ! area of v cells |
|---|
| 34 | !$OMP THREADPRIVATE(airev) |
|---|
| 35 | |
|---|
| 36 | ! Local variables for horiz mass flux in slab |
|---|
| 37 | LOGICAL,SAVE :: alpha_var |
|---|
| 38 | !$OMP THREADPRIVATE(alpha_var) |
|---|
| 39 | LOGICAL,SAVE :: slab_upstream |
|---|
| 40 | !$OMP THREADPRIVATE(slab_upstream) |
|---|
| 41 | |
|---|
| 42 | REAL,SAVE,ALLOCATABLE :: zmasqu(:) ! continent mask for zonal mass flux |
|---|
| 43 | !$OMP THREADPRIVATE(zmasqu) |
|---|
| 44 | REAL,SAVE,ALLOCATABLE :: zmasqv(:) ! continent mask for merid mass flux |
|---|
| 45 | !$OMP THREADPRIVATE(zmasqv) |
|---|
| 46 | REAL,SAVE,ALLOCATABLE :: unsfv(:) ! 1/f, v points |
|---|
| 47 | !$OMP THREADPRIVATE(unsfv) |
|---|
| 48 | REAL,SAVE,ALLOCATABLE :: unsev(:) ! 1/epsilon (drag) |
|---|
| 49 | !$OMP THREADPRIVATE(unsev) |
|---|
| 50 | REAL,SAVE,ALLOCATABLE :: unsfu(:) ! 1/F, u points |
|---|
| 51 | !$OMP THREADPRIVATE(unsfu) |
|---|
| 52 | REAL,SAVE,ALLOCATABLE :: unseu(:) |
|---|
| 53 | !$OMP THREADPRIVATE(unseu) |
|---|
| 54 | |
|---|
| 55 | ! Routines from dyn3d, valid on global dynamics grid only: |
|---|
| 56 | PRIVATE :: gr_fi_dyn, gr_dyn_fi ! to go between 1D nd 2D horiz grid |
|---|
| 57 | PRIVATE :: gr_scal_v,gr_v_scal,gr_scal_u ! change on 2D grid U,V, T points |
|---|
| 58 | PRIVATE :: grad,diverg |
|---|
| 59 | |
|---|
| 60 | CONTAINS |
|---|
| 61 | |
|---|
| 62 | SUBROUTINE ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez_,fext_,unsaire_,& |
|---|
| 63 | cu_,cuvsurcv_,cv_,cvusurcu_, & |
|---|
| 64 | aire_,apoln_,apols_, & |
|---|
| 65 | aireu_,airev_) |
|---|
| 66 | ! number of points in lon, lat |
|---|
| 67 | IMPLICIT NONE |
|---|
| 68 | ! Routine copies some geometry variables from the dynamical core |
|---|
| 69 | ! see global vars for meaning |
|---|
| 70 | INTEGER,INTENT(IN) :: ip1jm |
|---|
| 71 | INTEGER,INTENT(IN) :: ip1jmp1 |
|---|
| 72 | REAL,INTENT(IN) :: unsairez_(ip1jm) |
|---|
| 73 | REAL,INTENT(IN) :: fext_(ip1jm) |
|---|
| 74 | REAL,INTENT(IN) :: unsaire_(ip1jmp1) |
|---|
| 75 | REAL,INTENT(IN) :: cu_(ip1jmp1) |
|---|
| 76 | REAL,INTENT(IN) :: cuvsurcv_(ip1jm) |
|---|
| 77 | REAL,INTENT(IN) :: cv_(ip1jm) |
|---|
| 78 | REAL,INTENT(IN) :: cvusurcu_(ip1jmp1) |
|---|
| 79 | REAL,INTENT(IN) :: aire_(ip1jmp1) |
|---|
| 80 | REAL,INTENT(IN) :: apoln_ |
|---|
| 81 | REAL,INTENT(IN) :: apols_ |
|---|
| 82 | REAL,INTENT(IN) :: aireu_(ip1jmp1) |
|---|
| 83 | REAL,INTENT(IN) :: airev_(ip1jm) |
|---|
| 84 | |
|---|
| 85 | |
|---|
| 86 | ! Sanity check on dimensions |
|---|
| 87 | if ((ip1jm.ne.((nbp_lon+1)*(nbp_lat-1))).or. & |
|---|
| 88 | (ip1jmp1.ne.((nbp_lon+1)*nbp_lat))) then |
|---|
| 89 | write(*,*) "ini_slab_transp_geom Error: wrong array sizes" |
|---|
| 90 | stop |
|---|
| 91 | endif |
|---|
| 92 | ! Allocations could be done only on master process/thread... |
|---|
| 93 | allocate(unsairez(ip1jm)) |
|---|
| 94 | unsairez(:)=unsairez_(:) |
|---|
| 95 | allocate(fext(ip1jm)) |
|---|
| 96 | fext(:)=fext_(:) |
|---|
| 97 | allocate(unsaire(ip1jmp1)) |
|---|
| 98 | unsaire(:)=unsaire_(:) |
|---|
| 99 | allocate(cu(ip1jmp1)) |
|---|
| 100 | cu(:)=cu_(:) |
|---|
| 101 | allocate(cuvsurcv(ip1jm)) |
|---|
| 102 | cuvsurcv(:)=cuvsurcv_(:) |
|---|
| 103 | allocate(cv(ip1jm)) |
|---|
| 104 | cv(:)=cv_(:) |
|---|
| 105 | allocate(cvusurcu(ip1jmp1)) |
|---|
| 106 | cvusurcu(:)=cvusurcu_(:) |
|---|
| 107 | allocate(aire(ip1jmp1)) |
|---|
| 108 | aire(:)=aire_(:) |
|---|
| 109 | apoln=apoln_ |
|---|
| 110 | apols=apols_ |
|---|
| 111 | allocate(aireu(ip1jmp1)) |
|---|
| 112 | aireu(:)=aireu_(:) |
|---|
| 113 | allocate(airev(ip1jm)) |
|---|
| 114 | airev(:)=airev_(:) |
|---|
| 115 | |
|---|
| 116 | END SUBROUTINE ini_slab_transp_geom |
|---|
| 117 | |
|---|
| 118 | SUBROUTINE ini_slab_transp(zmasq) |
|---|
| 119 | |
|---|
| 120 | ! USE ioipsl_getin_p_mod, only: getin_p |
|---|
| 121 | USE IOIPSL, ONLY : getin |
|---|
| 122 | IMPLICIT NONE |
|---|
| 123 | |
|---|
| 124 | REAL zmasq(klon_glo) ! ocean / continent mask, 1=continent |
|---|
| 125 | REAL zmasq_2d((nbp_lon+1)*nbp_lat) |
|---|
| 126 | REAL ff((nbp_lon+1)*(nbp_lat-1)) ! Coriolis parameter |
|---|
| 127 | REAL eps ! epsilon friction timescale (s-1) |
|---|
| 128 | INTEGER :: slab_ekman |
|---|
| 129 | INTEGER i |
|---|
| 130 | INTEGER :: iim,iip1,jjp1,ip1jm,ip1jmp1 |
|---|
| 131 | |
|---|
| 132 | ! Some definition for grid size |
|---|
| 133 | ip1jm=(nbp_lon+1)*(nbp_lat-1) |
|---|
| 134 | ip1jmp1=(nbp_lon+1)*nbp_lat |
|---|
| 135 | iim=nbp_lon |
|---|
| 136 | iip1=nbp_lon+1 |
|---|
| 137 | jjp1=nbp_lat |
|---|
| 138 | ip1jm=(nbp_lon+1)*(nbp_lat-1) |
|---|
| 139 | ip1jmp1=(nbp_lon+1)*nbp_lat |
|---|
| 140 | |
|---|
| 141 | ! Define ocean / continent mask (no flux into continent cell) |
|---|
| 142 | allocate(zmasqu(ip1jmp1)) |
|---|
| 143 | allocate(zmasqv(ip1jm)) |
|---|
| 144 | zmasqu=1. |
|---|
| 145 | zmasqv=1. |
|---|
| 146 | |
|---|
| 147 | ! convert mask to 2D grid |
|---|
| 148 | CALL gr_fi_dyn(1,iip1,jjp1,zmasq,zmasq_2d) |
|---|
| 149 | ! put flux mask to 0 at boundaries of continent cells |
|---|
| 150 | DO i=1,ip1jmp1-1 |
|---|
| 151 | IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+1).GT.1e-5) THEN |
|---|
| 152 | zmasqu(i)=0. |
|---|
| 153 | ENDIF |
|---|
| 154 | END DO |
|---|
| 155 | DO i=iip1,ip1jmp1,iip1 |
|---|
| 156 | zmasqu(i)=zmasqu(i-iim) |
|---|
| 157 | END DO |
|---|
| 158 | DO i=1,ip1jm |
|---|
| 159 | IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+iip1).GT.1e-5) THEN |
|---|
| 160 | zmasqv(i)=0. |
|---|
| 161 | END IF |
|---|
| 162 | END DO |
|---|
| 163 | slab_ekman=2 |
|---|
| 164 | CALL getin("slab_ekman",slab_ekman) |
|---|
| 165 | ! Coriolis and friction for Ekman transport |
|---|
| 166 | IF (slab_ekman.GT.0) THEN |
|---|
| 167 | allocate(unsfv(ip1jm)) |
|---|
| 168 | allocate(unsev(ip1jm)) |
|---|
| 169 | allocate(unsfu(ip1jmp1)) |
|---|
| 170 | allocate(unseu(ip1jmp1)) |
|---|
| 171 | |
|---|
| 172 | eps=1e-5 ! Drag |
|---|
| 173 | CALL getin('slab_eps',eps) |
|---|
| 174 | print *,'epsilon=',eps |
|---|
| 175 | ff=fext*unsairez ! Coriolis |
|---|
| 176 | ! coefs to convert tau_x, tau_y to Ekman mass fluxes |
|---|
| 177 | ! on 2D grid v points |
|---|
| 178 | DO i=1,ip1jm |
|---|
| 179 | unsev(i)=eps/(ff(i)*ff(i)+eps**2) |
|---|
| 180 | unsfv(i)=ff(i)/(ff(i)*ff(i)+eps**2) |
|---|
| 181 | ENDDO |
|---|
| 182 | ! compute values on 2D u grid |
|---|
| 183 | CALL gr_v_scal(1,unsfv,unsfu) |
|---|
| 184 | CALL gr_v_scal(1,unsev,unseu) |
|---|
| 185 | END IF |
|---|
| 186 | |
|---|
| 187 | ! Other options for Ekman transport |
|---|
| 188 | ! Alpha variable? |
|---|
| 189 | alpha_var=.FALSE. |
|---|
| 190 | CALL getin('slab_alphav',alpha_var) |
|---|
| 191 | print *,'alpha variable',alpha_var |
|---|
| 192 | ! centered ou upstream scheme for meridional transport |
|---|
| 193 | slab_upstream=.FALSE. |
|---|
| 194 | CALL getin('slab_upstream',slab_upstream) |
|---|
| 195 | print *,'upstream slab scheme',slab_upstream |
|---|
| 196 | |
|---|
| 197 | END SUBROUTINE ini_slab_transp |
|---|
| 198 | |
|---|
| 199 | SUBROUTINE divgrad_phy(nlevs,temp,delta) |
|---|
| 200 | ! Computes temperature tendency due to horizontal diffusion : |
|---|
| 201 | ! T Laplacian, later multiplied by diffusion coef and time-step |
|---|
| 202 | |
|---|
| 203 | IMPLICIT NONE |
|---|
| 204 | |
|---|
| 205 | INTEGER, INTENT(IN) :: nlevs ! nlevs : slab layers |
|---|
| 206 | REAL, INTENT(IN) :: temp(klon_glo,nlevs) ! slab temperature |
|---|
| 207 | REAL , INTENT(OUT) :: delta(klon_glo,nlevs) ! temp laplacian (heat flux div.) |
|---|
| 208 | REAL :: delta_2d((nbp_lon+1)*nbp_lat,nlevs) |
|---|
| 209 | REAL ghx((nbp_lon+1)*nbp_lat,nlevs), ghy((nbp_lon+1)*(nbp_lat-1),nlevs) |
|---|
| 210 | INTEGER :: ll,iip1,jjp1 |
|---|
| 211 | |
|---|
| 212 | iip1=nbp_lon+1 |
|---|
| 213 | jjp1=nbp_lat |
|---|
| 214 | |
|---|
| 215 | ! transpose temp to 2D horiz. grid |
|---|
| 216 | CALL gr_fi_dyn(nlevs,iip1,jjp1,temp,delta_2d) |
|---|
| 217 | ! computes gradient (proportional to heat flx) |
|---|
| 218 | CALL grad(nlevs,delta_2d,ghx,ghy) |
|---|
| 219 | ! put flux to 0 at ocean / continent boundary |
|---|
| 220 | DO ll=1,nlevs |
|---|
| 221 | ghx(:,ll)=ghx(:,ll)*zmasqu |
|---|
| 222 | ghy(:,ll)=ghy(:,ll)*zmasqv |
|---|
| 223 | END DO |
|---|
| 224 | ! flux divergence |
|---|
| 225 | CALL diverg(nlevs,ghx,ghy,delta_2d) |
|---|
| 226 | ! laplacian back to 1D grid |
|---|
| 227 | CALL gr_dyn_fi(nlevs,iip1,jjp1,delta_2d,delta) |
|---|
| 228 | |
|---|
| 229 | RETURN |
|---|
| 230 | END SUBROUTINE divgrad_phy |
|---|
| 231 | |
|---|
| 232 | SUBROUTINE slab_ekman1(tx_phy,ty_phy,ts_phy,dt_phy) |
|---|
| 233 | ! 1.5 Layer Ekman transport temperature tendency |
|---|
| 234 | ! note : tendency dt later multiplied by (delta t)/(rho.H) |
|---|
| 235 | ! to convert from divergence of heat fluxes to T |
|---|
| 236 | |
|---|
| 237 | IMPLICIT NONE |
|---|
| 238 | |
|---|
| 239 | ! tx, ty : wind stress (different grids) |
|---|
| 240 | ! fluxm, fluz : mass *or heat* fluxes |
|---|
| 241 | ! dt : temperature tendency |
|---|
| 242 | INTEGER ij |
|---|
| 243 | |
|---|
| 244 | ! ts surface temp, td deep temp (diagnosed) |
|---|
| 245 | REAL ts_phy(klon_glo) |
|---|
| 246 | REAL ts((nbp_lon+1)*nbp_lat),td((nbp_lon+1)*nbp_lat) |
|---|
| 247 | ! zonal and meridional wind stress |
|---|
| 248 | REAL tx_phy(klon_glo),ty_phy(klon_glo) |
|---|
| 249 | REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat) |
|---|
| 250 | REAL txv((nbp_lon+1)*(nbp_lat-1)),tyv((nbp_lon+1)*(nbp_lat-1)) |
|---|
| 251 | ! zonal and meridional Ekman mass fluxes at u, v points (2D grid) |
|---|
| 252 | REAL fluxz((nbp_lon+1)*nbp_lat),fluxm((nbp_lon+1)*(nbp_lat-1)) |
|---|
| 253 | ! vertical and absolute mass fluxes (to estimate alpha) |
|---|
| 254 | REAL fluxv((nbp_lon+1)*nbp_lat),fluxt((nbp_lon+1)*(nbp_lat-1)) |
|---|
| 255 | ! temperature tendency |
|---|
| 256 | REAL dt((nbp_lon+1)*nbp_lat),dt_phy(klon_glo) |
|---|
| 257 | REAL alpha((nbp_lon+1)*nbp_lat) ! deep temperature coef |
|---|
| 258 | |
|---|
| 259 | INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1 |
|---|
| 260 | |
|---|
| 261 | ! Grid definitions |
|---|
| 262 | iim=nbp_lon |
|---|
| 263 | iip1=nbp_lon+1 |
|---|
| 264 | iip2=nbp_lon+2 |
|---|
| 265 | jjp1=nbp_lat |
|---|
| 266 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
|---|
| 267 | ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1 |
|---|
| 268 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
|---|
| 269 | |
|---|
| 270 | ! Convert taux,y to 2D scalar grid |
|---|
| 271 | ! Note: 2D grid size = iim*jjm. iip1=iim+1 |
|---|
| 272 | ! First and last points in zonal direction are the same |
|---|
| 273 | ! we use 1 index ij from 1 to (iim+1)*(jjm+1) |
|---|
| 274 | CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu) |
|---|
| 275 | CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu) |
|---|
| 276 | ! convert to u,v grid (Arakawa C) |
|---|
| 277 | ! Multiply by f or eps to get mass flux |
|---|
| 278 | ! Meridional mass flux |
|---|
| 279 | CALL gr_scal_v(1,txu,txv) ! wind stress at v points |
|---|
| 280 | CALL gr_scal_v(1,tyu,tyv) |
|---|
| 281 | fluxm=tyv*unsev-txv*unsfv ! in kg.s-1.m-1 (zonal distance) |
|---|
| 282 | ! Zonal mass flux |
|---|
| 283 | CALL gr_scal_u(1,txu,txu) ! wind stress at u points |
|---|
| 284 | CALL gr_scal_u(1,tyu,tyu) |
|---|
| 285 | fluxz=tyu*unsfu+txu*unseu |
|---|
| 286 | |
|---|
| 287 | ! Correct flux: continent mask and horiz grid size |
|---|
| 288 | ! multiply m-flux by mask and dx: flux in kg.s-1 |
|---|
| 289 | fluxm=fluxm*cv*cuvsurcv*zmasqv |
|---|
| 290 | ! multiply z-flux by mask and dy: flux in kg.s-1 |
|---|
| 291 | fluxz=fluxz*cu*cvusurcu*zmasqu |
|---|
| 292 | |
|---|
| 293 | ! Compute vertical and absolute mass flux (for variable alpha) |
|---|
| 294 | IF (alpha_var) THEN |
|---|
| 295 | DO ij=iip2,ip1jm |
|---|
| 296 | fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1) |
|---|
| 297 | fluxt(ij)=ABS(fluxz(ij))+ABS(fluxz(ij-1)) & |
|---|
| 298 | +ABS(fluxm(ij))+ABS(fluxm(ij-iip1)) |
|---|
| 299 | ENDDO |
|---|
| 300 | DO ij=iip1,ip1jmi1,iip1 |
|---|
| 301 | fluxt(ij+1)=fluxt(ij+iip1) |
|---|
| 302 | fluxv(ij+1)=fluxv(ij+iip1) |
|---|
| 303 | END DO |
|---|
| 304 | fluxt(1)=SUM(ABS(fluxm(1:iim))) |
|---|
| 305 | fluxt(ip1jmp1)=SUM(ABS(fluxm(ip1jm-iim:ip1jm-1))) |
|---|
| 306 | fluxv(1)=-SUM(fluxm(1:iim)) |
|---|
| 307 | fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1)) |
|---|
| 308 | fluxt=MAX(fluxt,1.e-10) |
|---|
| 309 | ENDIF |
|---|
| 310 | |
|---|
| 311 | ! Compute alpha coefficient. |
|---|
| 312 | ! Tdeep = Tsurf * alpha + 271.35 * (1-alpha) |
|---|
| 313 | IF (alpha_var) THEN |
|---|
| 314 | ! increase alpha (and Tdeep) in downwelling regions |
|---|
| 315 | ! and decrease in upwelling regions |
|---|
| 316 | ! to avoid "hot spots" where there is surface convergence |
|---|
| 317 | DO ij=iip2,ip1jm |
|---|
| 318 | alpha(ij)=2./3.-fluxv(ij)/fluxt(ij)/3. |
|---|
| 319 | ENDDO |
|---|
| 320 | alpha(1:iip1)=2./3.-fluxv(1)/fluxt(1)/3. |
|---|
| 321 | alpha(ip1jm+1:ip1jmp1)=2./3.-fluxv(ip1jmp1)/fluxt(ip1jmp1)/3. |
|---|
| 322 | ELSE |
|---|
| 323 | alpha(:)=2./3. |
|---|
| 324 | ! Tsurf-Tdeep ~ 10° in the Tropics |
|---|
| 325 | ENDIF |
|---|
| 326 | |
|---|
| 327 | ! Estimate deep temperature |
|---|
| 328 | CALL gr_fi_dyn(1,iip1,jjp1,ts_phy,ts) |
|---|
| 329 | DO ij=1,ip1jmp1 |
|---|
| 330 | td(ij)=271.35+(ts(ij)-271.35)*alpha(ij) |
|---|
| 331 | td(ij)=MIN(td(ij),ts(ij)) |
|---|
| 332 | END DO |
|---|
| 333 | |
|---|
| 334 | ! Meridional heat flux: multiply mass flux by (ts-td) |
|---|
| 335 | ! flux in K.kg.s-1 |
|---|
| 336 | IF (slab_upstream) THEN |
|---|
| 337 | ! upstream scheme to avoid hot spots |
|---|
| 338 | DO ij=1,ip1jm |
|---|
| 339 | IF (fluxm(ij).GE.0.) THEN |
|---|
| 340 | fluxm(ij)=fluxm(ij)*(ts(ij+iip1)-td(ij)) |
|---|
| 341 | ELSE |
|---|
| 342 | fluxm(ij)=fluxm(ij)*(ts(ij)-td(ij+iip1)) |
|---|
| 343 | END IF |
|---|
| 344 | END DO |
|---|
| 345 | ELSE |
|---|
| 346 | ! centered scheme better in mid-latitudes |
|---|
| 347 | DO ij=1,ip1jm |
|---|
| 348 | fluxm(ij)=fluxm(ij)*(ts(ij+iip1)+ts(ij)-td(ij)-td(ij+iip1))/2. |
|---|
| 349 | END DO |
|---|
| 350 | ENDIF |
|---|
| 351 | |
|---|
| 352 | ! Zonal heat flux |
|---|
| 353 | ! upstream scheme |
|---|
| 354 | DO ij=iip2,ip1jm |
|---|
| 355 | fluxz(ij)=fluxz(ij)*(ts(ij)+ts(ij+1)-td(ij+1)-td(ij))/2. |
|---|
| 356 | END DO |
|---|
| 357 | DO ij=iip1*2,ip1jmp1,iip1 |
|---|
| 358 | fluxz(ij)=fluxz(ij-iim) |
|---|
| 359 | END DO |
|---|
| 360 | |
|---|
| 361 | ! temperature tendency = divergence of heat fluxes |
|---|
| 362 | ! dt in K.s-1.kg.m-2 (T trend times mass/horiz surface) |
|---|
| 363 | DO ij=iip2,ip1jm |
|---|
| 364 | dt(ij)=(fluxz(ij-1)-fluxz(ij)+fluxm(ij)-fluxm(ij-iip1)) & |
|---|
| 365 | /aire(ij) ! aire : grid area |
|---|
| 366 | END DO |
|---|
| 367 | DO ij=iip1,ip1jmi1,iip1 |
|---|
| 368 | dt(ij+1)=dt(ij+iip1) |
|---|
| 369 | END DO |
|---|
| 370 | ! special treatment at the Poles |
|---|
| 371 | dt(1)=SUM(fluxm(1:iim))/apoln |
|---|
| 372 | dt(ip1jmp1)=-SUM(fluxm(ip1jm-iim:ip1jm-1))/apols |
|---|
| 373 | dt(2:iip1)=dt(1) |
|---|
| 374 | dt(ip1jm+1:ip1jmp1)=dt(ip1jmp1) |
|---|
| 375 | |
|---|
| 376 | ! tendencies back to 1D grid |
|---|
| 377 | CALL gr_dyn_fi(1,iip1,jjp1,dt,dt_phy) |
|---|
| 378 | |
|---|
| 379 | RETURN |
|---|
| 380 | END SUBROUTINE slab_ekman1 |
|---|
| 381 | |
|---|
| 382 | SUBROUTINE slab_ekman2(tx_phy,ty_phy,ts_phy,dt_phy) |
|---|
| 383 | ! Temperature tendency for 2-layers slab ocean |
|---|
| 384 | ! note : tendency dt later multiplied by (delta time)/(rho.H) |
|---|
| 385 | ! to convert from divergence of heat fluxes to T |
|---|
| 386 | |
|---|
| 387 | IMPLICIT NONE |
|---|
| 388 | |
|---|
| 389 | ! Here, temperature and flux variables are on 2 layers |
|---|
| 390 | INTEGER ij |
|---|
| 391 | |
|---|
| 392 | REAL tx_phy(klon_glo),ty_phy(klon_glo) |
|---|
| 393 | REAL txv((nbp_lon+1)*(nbp_lat-1)), tyv((nbp_lon+1)*(nbp_lat-1)) |
|---|
| 394 | REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat) |
|---|
| 395 | ! slab temperature on 1D, 2D grid |
|---|
| 396 | REAL ts_phy(klon_glo,2), ts((nbp_lon+1)*nbp_lat,2) |
|---|
| 397 | ! zonal and meridional Ekman mass flux at u, v points (2D grid) |
|---|
| 398 | REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1)) |
|---|
| 399 | ! vertical mass flux between the 2 layers |
|---|
| 400 | REAL fluxv((nbp_lon+1)*nbp_lat) |
|---|
| 401 | ! zonal and meridional heat fluxes |
|---|
| 402 | REAL fluxtz((nbp_lon+1)*nbp_lat,2) |
|---|
| 403 | REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2) |
|---|
| 404 | ! temperature tendency (in K.s-1.kg.m-2) |
|---|
| 405 | REAL dt((nbp_lon+1)*nbp_lat,2), dt_phy(klon_glo,2) |
|---|
| 406 | |
|---|
| 407 | INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1 |
|---|
| 408 | |
|---|
| 409 | ! Grid definitions |
|---|
| 410 | iim=nbp_lon |
|---|
| 411 | iip1=nbp_lon+1 |
|---|
| 412 | iip2=nbp_lon+2 |
|---|
| 413 | jjp1=nbp_lat |
|---|
| 414 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
|---|
| 415 | ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1 |
|---|
| 416 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
|---|
| 417 | |
|---|
| 418 | ! Convert taux,y to 2D scalar grid |
|---|
| 419 | CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu) |
|---|
| 420 | CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu) |
|---|
| 421 | ! Multiply taux,y by f or eps, and convert to 2D u,v grids |
|---|
| 422 | ! (Arakawa C grid) |
|---|
| 423 | ! Meridional flux |
|---|
| 424 | CALL gr_scal_v(1,txu,txv) ! wind stress at v points |
|---|
| 425 | CALL gr_scal_v(1,tyu,tyv) |
|---|
| 426 | fluxm=tyv*unsev-txv*unsfv |
|---|
| 427 | ! Zonal flux |
|---|
| 428 | CALL gr_scal_u(1,txu,txu) ! wind stress at u points |
|---|
| 429 | CALL gr_scal_u(1,tyu,tyu) |
|---|
| 430 | fluxz=tyu*unsfu+txu*unseu |
|---|
| 431 | |
|---|
| 432 | ! Convert temperature to 2D grid |
|---|
| 433 | CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts) |
|---|
| 434 | |
|---|
| 435 | ! Mass fluxes (apply continent mask, multiply by horiz grid dimension) |
|---|
| 436 | fluxm=fluxm*cv*cuvsurcv*zmasqv |
|---|
| 437 | fluxz=fluxz*cu*cvusurcu*zmasqu |
|---|
| 438 | ! Vertical mass flux from mass budget (divergence of horiz fluxes) |
|---|
| 439 | DO ij=iip2,ip1jm |
|---|
| 440 | fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1) |
|---|
| 441 | ENDDO |
|---|
| 442 | DO ij=iip1,ip1jmi1,iip1 |
|---|
| 443 | fluxv(ij+1)=fluxv(ij+iip1) |
|---|
| 444 | END DO |
|---|
| 445 | ! vertical mass flux at Poles |
|---|
| 446 | fluxv(1)=-SUM(fluxm(1:iim)) |
|---|
| 447 | fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1)) |
|---|
| 448 | fluxv=fluxv |
|---|
| 449 | |
|---|
| 450 | ! Meridional heat fluxes |
|---|
| 451 | DO ij=1,ip1jm |
|---|
| 452 | ! centered scheme |
|---|
| 453 | fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2. |
|---|
| 454 | fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2. |
|---|
| 455 | END DO |
|---|
| 456 | |
|---|
| 457 | ! Zonal heat fluxes |
|---|
| 458 | ! Schema upstream |
|---|
| 459 | DO ij=iip2,ip1jm |
|---|
| 460 | IF (fluxz(ij).GE.0.) THEN |
|---|
| 461 | fluxtz(ij,1)=fluxz(ij)*ts(ij,1) |
|---|
| 462 | fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2) |
|---|
| 463 | ELSE |
|---|
| 464 | fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1) |
|---|
| 465 | fluxtz(ij,2)=-fluxz(ij)*ts(ij,2) |
|---|
| 466 | ENDIF |
|---|
| 467 | ENDDO |
|---|
| 468 | DO ij=iip1*2,ip1jmp1,iip1 |
|---|
| 469 | fluxtz(ij,:)=fluxtz(ij-iim,:) |
|---|
| 470 | END DO |
|---|
| 471 | |
|---|
| 472 | ! Temperature tendency : |
|---|
| 473 | DO ij=iip2,ip1jm |
|---|
| 474 | ! divergence of horizontal heat fluxes |
|---|
| 475 | dt(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) & |
|---|
| 476 | +fluxtm(ij,:)-fluxtm(ij-iip1,:) |
|---|
| 477 | ! + vertical heat flux (mass flux * T, upstream scheme) |
|---|
| 478 | IF (fluxv(ij).GT.0.) THEN |
|---|
| 479 | dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,2) |
|---|
| 480 | dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,2) |
|---|
| 481 | ELSE |
|---|
| 482 | dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,1) |
|---|
| 483 | dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,1) |
|---|
| 484 | ENDIF |
|---|
| 485 | ! divide by cell area |
|---|
| 486 | dt(ij,:)=dt(ij,:)/aire(ij) |
|---|
| 487 | END DO |
|---|
| 488 | DO ij=iip1,ip1jmi1,iip1 |
|---|
| 489 | dt(ij+1,:)=dt(ij+iip1,:) |
|---|
| 490 | END DO |
|---|
| 491 | ! Pôles |
|---|
| 492 | dt(1,:)=SUM(fluxtm(1:iim,:),dim=1) |
|---|
| 493 | IF (fluxv(1).GT.0.) THEN |
|---|
| 494 | dt(1,1)=dt(1,1)+fluxv(1)*ts(1,2) |
|---|
| 495 | dt(1,2)=dt(1,2)-fluxv(1)*ts(1,2) |
|---|
| 496 | ELSE |
|---|
| 497 | dt(1,1)=dt(1,1)+fluxv(1)*ts(1,1) |
|---|
| 498 | dt(1,2)=dt(1,2)-fluxv(1)*ts(1,1) |
|---|
| 499 | ENDIF |
|---|
| 500 | dt(1,:)=dt(1,:)/apoln |
|---|
| 501 | dt(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1) |
|---|
| 502 | IF (fluxv(ip1jmp1).GT.0.) THEN |
|---|
| 503 | dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,2) |
|---|
| 504 | dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,2) |
|---|
| 505 | ELSE |
|---|
| 506 | dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,1) |
|---|
| 507 | dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,1) |
|---|
| 508 | ENDIF |
|---|
| 509 | dt(ip1jmp1,:)=dt(ip1jmp1,:)/apols |
|---|
| 510 | dt(2:iip1,1)=dt(1,1) |
|---|
| 511 | dt(2:iip1,2)=dt(1,2) |
|---|
| 512 | dt(ip1jm+1:ip1jmp1,1)=dt(ip1jmp1,1) |
|---|
| 513 | dt(ip1jm+1:ip1jmp1,2)=dt(ip1jmp1,2) |
|---|
| 514 | |
|---|
| 515 | ! T tendency back to 1D grid... |
|---|
| 516 | CALL gr_dyn_fi(2,iip1,jjp1,dt,dt_phy) |
|---|
| 517 | |
|---|
| 518 | RETURN |
|---|
| 519 | END SUBROUTINE slab_ekman2 |
|---|
| 520 | |
|---|
| 521 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 522 | |
|---|
| 523 | SUBROUTINE gr_fi_dyn(nfield,im,jm,pfi,pdyn) |
|---|
| 524 | ! Transfer a variable from 1D "physics" grid to 2D "dynamics" grid |
|---|
| 525 | IMPLICIT NONE |
|---|
| 526 | |
|---|
| 527 | INTEGER,INTENT(IN) :: im,jm,nfield |
|---|
| 528 | REAL,INTENT(IN) :: pfi(klon_glo,nfield) ! on 1D grid |
|---|
| 529 | REAL,INTENT(OUT) :: pdyn(im,jm,nfield) ! on 2D grid |
|---|
| 530 | |
|---|
| 531 | INTEGER :: i,j,ifield,ig |
|---|
| 532 | |
|---|
| 533 | DO ifield=1,nfield |
|---|
| 534 | ! Handle poles |
|---|
| 535 | DO i=1,im |
|---|
| 536 | pdyn(i,1,ifield)=pfi(1,ifield) |
|---|
| 537 | pdyn(i,jm,ifield)=pfi(klon_glo,ifield) |
|---|
| 538 | ENDDO |
|---|
| 539 | ! Other points |
|---|
| 540 | DO j=2,jm-1 |
|---|
| 541 | ig=2+(j-2)*(im-1) |
|---|
| 542 | CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1) |
|---|
| 543 | pdyn(im,j,ifield)=pdyn(1,j,ifield) |
|---|
| 544 | ENDDO |
|---|
| 545 | ENDDO ! of DO ifield=1,nfield |
|---|
| 546 | |
|---|
| 547 | END SUBROUTINE gr_fi_dyn |
|---|
| 548 | |
|---|
| 549 | SUBROUTINE gr_dyn_fi(nfield,im,jm,pdyn,pfi) |
|---|
| 550 | ! Transfer a variable from 2D "dynamics" grid to 1D "physics" grid |
|---|
| 551 | IMPLICIT NONE |
|---|
| 552 | |
|---|
| 553 | INTEGER,INTENT(IN) :: im,jm,nfield |
|---|
| 554 | REAL,INTENT(IN) :: pdyn(im,jm,nfield) ! on 2D grid |
|---|
| 555 | REAL,INTENT(OUT) :: pfi(klon_glo,nfield) ! on 1D grid |
|---|
| 556 | |
|---|
| 557 | INTEGER j,ifield,ig |
|---|
| 558 | |
|---|
| 559 | ! Sanity check: |
|---|
| 560 | IF(klon_glo.NE.2+(jm-2)*(im-1)) THEN |
|---|
| 561 | WRITE(*,*) "gr_dyn_fi error, wrong sizes" |
|---|
| 562 | STOP |
|---|
| 563 | ENDIF |
|---|
| 564 | |
|---|
| 565 | ! Handle poles |
|---|
| 566 | CALL SCOPY(nfield,pdyn,im*jm,pfi,klon_glo) |
|---|
| 567 | CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(klon_glo,1),klon_glo) |
|---|
| 568 | ! Other points |
|---|
| 569 | DO ifield=1,nfield |
|---|
| 570 | DO j=2,jm-1 |
|---|
| 571 | ig=2+(j-2)*(im-1) |
|---|
| 572 | CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1) |
|---|
| 573 | ENDDO |
|---|
| 574 | ENDDO |
|---|
| 575 | |
|---|
| 576 | END SUBROUTINE gr_dyn_fi |
|---|
| 577 | |
|---|
| 578 | SUBROUTINE grad(klevel,pg,pgx,pgy) |
|---|
| 579 | ! compute the covariant components pgx,pgy of the gradient of pg |
|---|
| 580 | ! pgx = d(pg)/dx * delta(x) = delta(pg) |
|---|
| 581 | IMPLICIT NONE |
|---|
| 582 | |
|---|
| 583 | INTEGER,INTENT(IN) :: klevel |
|---|
| 584 | REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel) |
|---|
| 585 | REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel) |
|---|
| 586 | REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel) |
|---|
| 587 | |
|---|
| 588 | INTEGER :: l,ij |
|---|
| 589 | INTEGER :: iim,iip1,ip1jm,ip1jmp1 |
|---|
| 590 | |
|---|
| 591 | iim=nbp_lon |
|---|
| 592 | iip1=nbp_lon+1 |
|---|
| 593 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
|---|
| 594 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
|---|
| 595 | |
|---|
| 596 | DO l=1,klevel |
|---|
| 597 | DO ij=1,ip1jmp1-1 |
|---|
| 598 | pgx(ij,l)=pg(ij+1,l)-pg(ij,l) |
|---|
| 599 | ENDDO |
|---|
| 600 | ! correction for pgx(ip1,j,l) ... |
|---|
| 601 | ! ... pgx(iip1,j,l)=pgx(1,j,l) ... |
|---|
| 602 | DO ij=iip1,ip1jmp1,iip1 |
|---|
| 603 | pgx(ij,l)=pgx(ij-iim,l) |
|---|
| 604 | ENDDO |
|---|
| 605 | DO ij=1,ip1jm |
|---|
| 606 | pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l) |
|---|
| 607 | ENDDO |
|---|
| 608 | ENDDO |
|---|
| 609 | |
|---|
| 610 | END SUBROUTINE grad |
|---|
| 611 | |
|---|
| 612 | SUBROUTINE diverg(klevel,x,y,div) |
|---|
| 613 | ! computes the divergence of a vector field of components |
|---|
| 614 | ! x,y. x and y being covariant components |
|---|
| 615 | IMPLICIT NONE |
|---|
| 616 | |
|---|
| 617 | INTEGER,INTENT(IN) :: klevel |
|---|
| 618 | REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel) |
|---|
| 619 | REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel) |
|---|
| 620 | REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel) |
|---|
| 621 | |
|---|
| 622 | INTEGER :: l,ij |
|---|
| 623 | INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1 |
|---|
| 624 | |
|---|
| 625 | REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1) |
|---|
| 626 | REAL :: sumypn,sumyps |
|---|
| 627 | REAL,EXTERNAL :: SSUM |
|---|
| 628 | |
|---|
| 629 | iim=nbp_lon |
|---|
| 630 | iip1=nbp_lon+1 |
|---|
| 631 | iip2=nbp_lon+2 |
|---|
| 632 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
|---|
| 633 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
|---|
| 634 | ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1 |
|---|
| 635 | |
|---|
| 636 | DO l=1,klevel |
|---|
| 637 | DO ij=iip2,ip1jm-1 |
|---|
| 638 | div(ij+1,l)= & |
|---|
| 639 | cvusurcu(ij+1)*x(ij+1,l)-cvusurcu(ij)*x(ij,l)+ & |
|---|
| 640 | cuvsurcv(ij-iim)*y(ij-iim,l)-cuvsurcv(ij+1)*y(ij+1,l) |
|---|
| 641 | ENDDO |
|---|
| 642 | ! correction for div(1,j,l) ... |
|---|
| 643 | ! ... div(1,j,l)= div(iip1,j,l) ... |
|---|
| 644 | DO ij=iip2,ip1jm,iip1 |
|---|
| 645 | div(ij,l)=div(ij+iim,l) |
|---|
| 646 | ENDDO |
|---|
| 647 | ! at the poles |
|---|
| 648 | DO ij=1,iim |
|---|
| 649 | aiy1(ij)=cuvsurcv(ij)*y(ij,l) |
|---|
| 650 | aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l) |
|---|
| 651 | ENDDO |
|---|
| 652 | sumypn=SSUM(iim,aiy1,1)/apoln |
|---|
| 653 | sumyps=SSUM(iim,aiy2,1)/apols |
|---|
| 654 | DO ij=1,iip1 |
|---|
| 655 | div(ij,l)=-sumypn |
|---|
| 656 | div(ij+ip1jm,l)=sumyps |
|---|
| 657 | ENDDO |
|---|
| 658 | ! End (poles) |
|---|
| 659 | ENDDO ! of DO l=1,klevel |
|---|
| 660 | |
|---|
| 661 | !!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 ) |
|---|
| 662 | DO l=1,klevel |
|---|
| 663 | DO ij=iip2,ip1jm |
|---|
| 664 | div(ij,l)=div(ij,l)*unsaire(ij) |
|---|
| 665 | ENDDO |
|---|
| 666 | ENDDO |
|---|
| 667 | |
|---|
| 668 | END SUBROUTINE diverg |
|---|
| 669 | |
|---|
| 670 | SUBROUTINE gr_v_scal(nx,x_v,x_scal) |
|---|
| 671 | ! convert values from v points to scalar points on C-grid |
|---|
| 672 | ! used to compute unsfu, unseu (u points, but depends only on latitude) |
|---|
| 673 | IMPLICIT NONE |
|---|
| 674 | |
|---|
| 675 | INTEGER,INTENT(IN) :: nx ! number of levels or fields |
|---|
| 676 | REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx) |
|---|
| 677 | REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx) |
|---|
| 678 | |
|---|
| 679 | INTEGER :: l,ij |
|---|
| 680 | INTEGER :: iip1,iip2,ip1jm,ip1jmp1 |
|---|
| 681 | |
|---|
| 682 | iip1=nbp_lon+1 |
|---|
| 683 | iip2=nbp_lon+2 |
|---|
| 684 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
|---|
| 685 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
|---|
| 686 | |
|---|
| 687 | DO l=1,nx |
|---|
| 688 | DO ij=iip2,ip1jm |
|---|
| 689 | x_scal(ij,l)= & |
|---|
| 690 | (airev(ij-iip1)*x_v(ij-iip1,l)+airev(ij)*x_v(ij,l)) & |
|---|
| 691 | /(airev(ij-iip1)+airev(ij)) |
|---|
| 692 | ENDDO |
|---|
| 693 | DO ij=1,iip1 |
|---|
| 694 | x_scal(ij,l)=0. |
|---|
| 695 | ENDDO |
|---|
| 696 | DO ij=ip1jm+1,ip1jmp1 |
|---|
| 697 | x_scal(ij,l)=0. |
|---|
| 698 | ENDDO |
|---|
| 699 | ENDDO |
|---|
| 700 | |
|---|
| 701 | END SUBROUTINE gr_v_scal |
|---|
| 702 | |
|---|
| 703 | SUBROUTINE gr_scal_v(nx,x_scal,x_v) |
|---|
| 704 | ! convert values from scalar points to v points on C-grid |
|---|
| 705 | ! used to compute wind stress at V points |
|---|
| 706 | IMPLICIT NONE |
|---|
| 707 | |
|---|
| 708 | INTEGER,INTENT(IN) :: nx ! number of levels or fields |
|---|
| 709 | REAL,INTENT(OUT) :: x_v((nbp_lon+1)*(nbp_lat-1),nx) |
|---|
| 710 | REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx) |
|---|
| 711 | |
|---|
| 712 | INTEGER :: l,ij |
|---|
| 713 | INTEGER :: iip1,ip1jm |
|---|
| 714 | |
|---|
| 715 | iip1=nbp_lon+1 |
|---|
| 716 | ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm |
|---|
| 717 | |
|---|
| 718 | DO l=1,nx |
|---|
| 719 | DO ij=1,ip1jm |
|---|
| 720 | x_v(ij,l)= & |
|---|
| 721 | (cu(ij)*cvusurcu(ij)*x_scal(ij,l)+ & |
|---|
| 722 | cu(ij+iip1)*cvusurcu(ij+iip1)*x_scal(ij+iip1,l)) & |
|---|
| 723 | /(cu(ij)*cvusurcu(ij)+cu(ij+iip1)*cvusurcu(ij+iip1)) |
|---|
| 724 | ENDDO |
|---|
| 725 | ENDDO |
|---|
| 726 | |
|---|
| 727 | END SUBROUTINE gr_scal_v |
|---|
| 728 | |
|---|
| 729 | SUBROUTINE gr_scal_u(nx,x_scal,x_u) |
|---|
| 730 | ! convert values from scalar points to U points on C-grid |
|---|
| 731 | ! used to compute wind stress at U points |
|---|
| 732 | IMPLICIT NONE |
|---|
| 733 | |
|---|
| 734 | INTEGER,INTENT(IN) :: nx |
|---|
| 735 | REAL,INTENT(OUT) :: x_u((nbp_lon+1)*nbp_lat,nx) |
|---|
| 736 | REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx) |
|---|
| 737 | |
|---|
| 738 | INTEGER :: l,ij |
|---|
| 739 | INTEGER :: iip1,jjp1,ip1jmp1 |
|---|
| 740 | |
|---|
| 741 | iip1=nbp_lon+1 |
|---|
| 742 | jjp1=nbp_lat |
|---|
| 743 | ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1 |
|---|
| 744 | |
|---|
| 745 | DO l=1,nx |
|---|
| 746 | DO ij=1,ip1jmp1-1 |
|---|
| 747 | x_u(ij,l)= & |
|---|
| 748 | (aire(ij)*x_scal(ij,l)+aire(ij+1)*x_scal(ij+1,l)) & |
|---|
| 749 | /(aire(ij)+aire(ij+1)) |
|---|
| 750 | ENDDO |
|---|
| 751 | ENDDO |
|---|
| 752 | |
|---|
| 753 | CALL SCOPY(nx*jjp1,x_u(1,1),iip1,x_u(iip1,1),iip1) |
|---|
| 754 | |
|---|
| 755 | END SUBROUTINE gr_scal_u |
|---|
| 756 | |
|---|
| 757 | END MODULE slab_heat_transp_mod |
|---|
