Changeset 3432
- Timestamp:
- 2012-07-11T13:22:58+02:00 (13 years ago)
- Location:
- branches/2011/DEV_r2739_STFC_dCSE
- Files:
-
- 17 added
- 1 deleted
- 45 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/CONFIG/GYRE/EXP00/namelist
r2715 r3432 27 27 cn_exp = "GYRE" ! experience name 28 28 nn_it000 = 1 ! first time step 29 nn_itend = 4320 ! last time step29 nn_itend = 10 ! 4320 ! last time step 30 30 nn_date0 = 010101 ! initial calendar date yymmdd (used if nn_rstctl=1) 31 31 nn_leapy = 30 ! Leap year calendar (1) or not (0) … … 679 679 &namsol ! elliptic solver / island / free surface 680 680 !----------------------------------------------------------------------- 681 nn_solv = 2! elliptic solver: =1 preconditioned conjugate gradient (pcg)681 nn_solv = 1 ! elliptic solver: =1 preconditioned conjugate gradient (pcg) 682 682 ! =2 successive-over-relaxation (sor) 683 683 nn_sol_arp = 0 ! absolute/relative (0/1) precision convergence test -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/CONFIG/GYRE/cpp_GYRE.fcm
r2670 r3432 1 bld::tool::fppkeys key_gyre key_dynspg_flt key_ldfslp key_zdftke key_vectopt_loop key_iomput1 bld::tool::fppkeys key_gyre key_dynspg_flt key_ldfslp key_zdftke key_iomput key_mpp_mpi key_mpp_rkpart -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
r3211 r3432 37 37 PUBLIC bdy_dta_alloc ! routine called by bdy_init.F90 38 38 39 INTEGER :: numbdyt, numbdyu, numbdyv 40 INTEGER :: ntimes_bdy 41 INTEGER :: nbdy_b, nbdy_a 42 INTEGER :: numbdyt_bt, numbdyu_bt, numbdyv_bt 43 INTEGER :: ntimes_bdy_bt 44 INTEGER :: nbdy_b_bt, nbdy_a_bt 45 46 INTEGER, DIMENSION (jpbtime) :: istep, istep_bt 47 48 REAL(wp) :: zoffset 39 INTEGER :: numbdyt, numbdyu, numbdyv ! logical units for T-, U-, & V-points data file, resp. 40 INTEGER :: ntimes_bdy ! exact number of time dumps in data files 41 INTEGER :: nbdy_b, nbdy_a ! record of bdy data file for before and after time step 42 INTEGER :: numbdyt_bt, numbdyu_bt, numbdyv_bt ! logical unit for T-, U- & V-points data file, resp. 43 INTEGER :: ntimes_bdy_bt ! exact number of time dumps in data files 44 INTEGER :: nbdy_b_bt, nbdy_a_bt ! record of bdy data file for before and after time step 45 46 INTEGER, DIMENSION (jpbtime) :: istep, istep_bt ! time array in seconds in each data file 47 48 REAL(wp) :: zoffset ! time offset between time origin in file & start time of model run 49 49 50 50 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tbdydta, sbdydta ! time interpolated values of T and S bdy data … … 118 118 REAL(wp) :: dayjul0, zdayjulini 119 119 REAL(wp), DIMENSION(jpbtime) :: zstepr ! REAL time array from data files 120 !! DCSE_NEMO: do not ftrans! Beware! 120 !! DCSE_NEMO: do not ftrans! Beware! ARP - why not? 121 !!!!! zdta :I :I :z 121 122 REAL(wp), DIMENSION(jpbdta,1,jpk) :: zdta ! temporary array for data fields 122 123 !!--------------------------------------------------------------------------- … … 215 216 zoffset = (dayjul0-zdayjulini)*86400 216 217 ! 217 7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 218 ! ARP - Cannot have 'string literals' in an input format with the Cray compiler 219 ! so replace with corresponding number of spaces instead. 220 !7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 221 7000 FORMAT(14x, I4.4,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2) 218 222 219 223 !! TO BE DONE... Check consistency between calendar from file … … 399 403 IF(lwp) WRITE(numout,*) 'Dim size for votemper is ', nblendta(igrd) 400 404 ipi = nblendta(igrd) 401 CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 405 ! DCSE_NEMO: Call iom_get and explicitly flag that array has 406 ! not been ftrans'd 407 CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', & 408 zdta(1:ipi,1:ipj,1:ipk), nbdy_a, lzfirst=.FALSE. ) 402 409 ! 403 410 DO ib = 1, nblen(igrd) … … 414 421 IF(lwp) WRITE(numout,*) 'Dim size for vosaline is ', nblendta(igrd) 415 422 ipi = nblendta(igrd) 416 CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 423 ! DCSE_NEMO: pass optional flag to signal that zdta is NOT 424 ! ftrans'd 425 CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', & 426 zdta(1:ipi,1:ipj,1:ipk), nbdy_a, lzfirst=.FALSE. ) 417 427 ! 418 428 DO ib = 1, nblen(igrd) … … 595 605 igrd = 1 ! temperature 596 606 ipi = nblendta(igrd) 597 CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 607 ! DCSE_NEMO: pass optional flag to signal that zdta is NOT 608 ! ftrans'd 609 CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', & 610 zdta(1:ipi,1:ipj,1:ipk), nbdy_a, lzfirst=.FALSE. ) 598 611 DO ib = 1, nblen(igrd) 599 612 DO ik = 1, jpkm1 … … 604 617 igrd = 1 ! salinity 605 618 ipi = nblendta(igrd) 606 CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 619 ! DCSE_NEMO: pass optional flag to signal that zdta is NOT 620 ! ftrans'd 621 CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', & 622 zdta(1:ipi,1:ipj,1:ipk), nbdy_a, lzfirst=.FALSE. ) 607 623 DO ib = 1, nblen(igrd) 608 624 DO ik = 1, jpkm1 … … 732 748 END SUBROUTINE bdy_dta_frs 733 749 750 !FTRANS CLEAR 734 751 735 752 SUBROUTINE bdy_dta_fla( kt, jit, icycl ) … … 872 889 zoffset = (dayjul0-zdayjulini)*86400 873 890 ! 874 875 7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 891 ! ARP - Cannot have 'string literals' in an input format with the Cray compiler 892 ! so replace with corresponding number of spaces instead. 893 !7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 894 7000 FORMAT(14x, I4.4,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2) 876 895 877 896 !! TO BE DONE... Check consistency between calendar from file -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90
r3211 r3432 26 26 USE bdytides ! for tidal harmonic forcing at boundary 27 27 USE in_out_manager ! 28 USE timing, ONLY: timing_start, timing_stop 28 29 29 30 IMPLICIT NONE … … 126 127 !! continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164. 127 128 !!---------------------------------------------------------------------- 129 USE exchmod, Only: add_exch, bound_exch_list 128 130 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh 129 131 … … 133 135 REAL(wp) :: zforc ! temporary scalar 134 136 !!---------------------------------------------------------------------- 137 138 CALL timing_start('bdy_dyn_fla') 135 139 136 140 ! ---------------------------------! … … 161 165 zcorr = - flagu(jb) * SQRT( grav * hur_e(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 162 166 zforc = ubtbdy(jb) + utide(jb) 163 ua_e(ii,ij) = zforc + zcorr * umask(ii,ij,1) 167 #if defined key_z_first 168 ua_e(ii,ij) = zforc + zcorr * umask_1(ii,ij) 169 #else 170 ua_e(ii,ij) = zforc + zcorr * umask(ii,ij,1) 171 #endif 164 172 END DO 165 173 ! … … 174 182 zcorr = - flagv(jb) * SQRT( grav * hvr_e(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 175 183 zforc = vbtbdy(jb) + vtide(jb) 184 #if defined key_z_first 185 va_e(ii,ij) = zforc + zcorr * vmask_1(ii,ij) 186 #else 176 187 va_e(ii,ij) = zforc + zcorr * vmask(ii,ij,1) 177 END DO 188 #endif 189 END DO 190 191 ! 192 #if defined key_mpp_rkpart 193 ! Is quicker to pack halo-swaps for 2D fields than to do them 194 ! individually 195 CALL add_exch(jpreci,'U',Iminus,Jminus,Iplus,Jplus,ua_e,isgn=-1) 196 CALL add_exch(jpreci,'V',Iminus,Jminus,Iplus,Jplus,va_e,isgn=-1) 197 CALL bound_exch_list() 198 #else 178 199 CALL lbc_lnk( ua_e, 'U', -1. ) ! Boundary points should be updated 179 200 CALL lbc_lnk( va_e, 'V', -1. ) ! 201 #endif 180 202 ! 181 203 ENDIF ! ln_dyn_fla .or. ln_tides 204 ! 205 CALL timing_stop('bdy_dyn_fla','section') 182 206 ! 183 207 END SUBROUTINE bdy_dyn_fla -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r3211 r3432 622 622 zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1) 623 623 624 IF(ln_ctl) CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' dw/nmln : ', tab2d_2=tmask_1, clinfo2=' dw/tm_1 : ', ovlap=1 ) 624 #if defined key_z_first 625 IF(ln_ctl) CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' dw/nmln : ', & 626 tab2d_2=tmask_1, clinfo2=' dw/tm_1 : ', & 627 ovlap=1 ) 628 #else 629 IF(ln_ctl) CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' dw/nmln : ', & 630 tab2d_2=tmask(:,:,1), clinfo2=' dw/tm_1 : ', & 631 ovlap=1 ) 632 #endif 625 633 IF(ln_ctl) CALL prt_ctl( tab2d_1=zw2d, clinfo1=' dw/zw2d : ', ovlap=1 ) 626 634 -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r3211 r3432 91 91 INTEGER, PUBLIC :: nbse, nbsw !: logical of south east & south west processor 92 92 INTEGER, PUBLIC :: nidom !: ??? 93 INTEGER, PUBLIC :: nwidthmax !: width of the widest northern domain inc. halos 93 94 94 95 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mig !: local ==> global domain i-index -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90
r3211 r3432 75 75 !! - mj0, , mj1 : 76 76 !!---------------------------------------------------------------------- 77 USE mapcomm_mod, Only: iesub, jesub 78 !! * Local declarations 77 79 INTEGER :: ji, jj ! dummy loop argument 78 80 !!---------------------------------------------------------------------- … … 81 83 ! ! Local domain ! 82 84 ! ! ============== ! 85 !!$#if defined key_mpp_rkpart 86 !!$ WRITE(*,*) 'ARPDBG: ',narea,': dom_glo: nld{i,j} = ',nldi,nldj 87 !!$ WRITE(*,*) 'ARPDBG: ',narea,': dom_glo: nlc{i,j} = ',nlci,nlcj 88 !!$ WRITE(*,*) 'ARPDBG: ',narea,': dom_glo: n{i,j}mpp = ',nimpp, njmpp 89 !!$ WRITE(*,*) 'ARPDBG: ',narea,': dom_glo: jp{i,j}zoom = ',jpizoom, jpjzoom 90 !!$ DO ji = 1, jpi, 1 91 !!$ mig(ji) = ji - nldi + jpizoom - 1 + nimpp - 1 92 !!$! WRITE(*,FMT="('ARPDBG, mig(',I3,') = ',I3)") & 93 !!$! ji, mig(ji) 94 !!$ END DO 95 !!$ DO jj = 1, jpj, 1 96 !!$ mjg(jj) = jj - nldj + jpjzoom - 1 + njmpp - 1 97 !!$! WRITE(*,FMT="('ARPDBG, mjg(',I3,') = ',I3)") & 98 !!$! jj, mjg(jj) 99 !!$ END DO 100 !!$#else 83 101 DO ji = 1, jpi ! local domain indices ==> data domain indices 84 102 mig(ji) = ji + jpizoom - 1 + nimpp - 1 … … 87 105 mjg(jj) = jj + jpjzoom - 1 + njmpp - 1 88 106 END DO 107 !!$#endif 89 108 ! 90 109 ! ! data domain indices ==> local domain indices 91 110 ! ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the 92 111 ! !local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain. 112 #if defined key_mpp_rkpart 113 mi0(1:nimpp-1) = 1 !nldi 114 DO ji = 0,iesub-1,1 115 mi0(nimpp+ji) = nldi + ji 116 END DO 117 ! Halo columns (if any) 118 IF(nimpp+iesub <= jpidta) mi0(nimpp+iesub) = jpi+1 119 IF(nimpp+iesub+1 <= jpidta) mi0(nimpp+iesub+1:)=jpi+1 120 121 mi1(1:nimpp-2) = 0 122 IF(nimpp > 1)mi1(nimpp - 1) = 0 123 DO ji = 0,iesub-1,1 124 mi1(nimpp+ji) = nldi + ji 125 END DO 126 ! Halo columns (if any) 127 IF(nimpp+iesub <= jpidta) mi1(nimpp+iesub) = jpi 128 IF(nimpp+iesub+1 <= jpidta) mi1(nimpp+iesub+1:)= jpi 129 130 mj0(:njmpp-1) = 1 ! nldj 131 DO jj = 0,jesub-1,1 132 mj0(njmpp+jj) = nldj + jj 133 END DO 134 ! Halo rows (if any) 135 IF(njmpp+jesub <= jpjdta) mj0(njmpp+jesub) = jpj+1 136 IF(njmpp+jesub+1 <= jpjdta) mj0(njmpp+jesub+1:)=jpj+1 137 138 mj1(1:njmpp-2) = 0 139 IF(njmpp > 1)mj1(njmpp - 1) = 0 140 DO jj = 0,jesub-1,1 141 mj1(njmpp+jj) = nldj + jj 142 END DO 143 ! Halo rows (if any) 144 IF(njmpp+jesub <= jpjdta) mj1(njmpp+jesub) = jpj 145 IF(njmpp+jesub+1 <= jpjdta) mj1(njmpp+jesub+1:)= jpj 146 #else 93 147 DO ji = 1, jpidta 94 148 mi0(ji) = MAX( 1, MIN( ji - jpizoom + 1 - nimpp + 1, jpi+1 ) ) … … 99 153 mj1(jj) = MAX( 0, MIN( jj - jpjzoom + 1 - njmpp + 1, jpj ) ) 100 154 END DO 155 #endif 101 156 102 157 IF(lwp) THEN ! control print -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r3211 r3432 156 156 ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 157 157 ! for ssh and scale factors 158 zs_t (:,:) = e1t(:,:) * e2t(:,:)159 zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)160 zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)158 zs_t (:,:) = e1t(:,:) * e2t(:,:) 159 zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) ) 160 zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) ) 161 161 162 162 DO jj = 1, jpjm1 ! initialise before and now Sea Surface Height at u-, v-, f-points -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r3226 r3432 54 54 ! ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 55 55 REAL(wp) :: rn_hc = 150._wp ! Critical depth for s-sigma coordinates 56 56 PUBLIC rn_sbot_min, rn_sbot_max, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb,rn_hc 57 57 !! * Control permutation of array indices 58 58 # include "oce_ftrans.h90" … … 82 82 !! ln_zco=T z-coordinate 83 83 !! ln_zps=T z-coordinate with partial steps 84 !! ln_ sco=T s-coordinate84 !! ln_zco=T s-coordinate 85 85 !! 86 86 !! ** Action : define gdep., e3., mbathy and bathy … … 395 395 mbathy(:,:) = 0 ! set to zero extra halo points 396 396 bathy (:,:) = 0._wp ! (require for mpp case) 397 #if defined key_mpp_rkpart 398 DO jj = nldj, nlcj ! interior values 399 DO ji = nldi, nlci 400 #else 397 401 DO jj = 1, nlcj ! interior values 398 402 DO ji = 1, nlci 403 #endif 399 404 mbathy(ji,jj) = idta( mig(ji), mjg(jj) ) 400 405 bathy (ji,jj) = zdta( mig(ji), mjg(jj) ) … … 593 598 !! - update bathy : meter bathymetry (in meters) 594 599 !!---------------------------------------------------------------------- 595 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 596 USE wrk_nemo, ONLY: zbathy => wrk_2d_1 600 USE mapcomm_mod, ONLY: trimmed, nidx,eidx,sidx,widx 601 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 602 USE wrk_nemo, ONLY: zbathy => wrk_2d_1 597 603 !! 598 604 INTEGER :: ji, jj, jl ! dummy loop indices … … 636 642 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed' 637 643 ENDIF 644 638 645 IF( lk_mpp ) THEN 639 646 zbathy(:,:) = FLOAT( mbathy(:,:) ) … … 646 653 IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio 647 654 IF( lk_mpp ) THEN 648 IF( nbondi == -1 .OR. nbondi == 2) THEN655 IF( (nbondi == -1 .OR. nbondi == 2) .AND. (.NOT. trimmed(widx,narea) ) ) THEN 649 656 IF( jperio /= 1 ) mbathy(1,:) = 0 650 657 ENDIF 651 IF( nbondi == 1 .OR. nbondi == 2) THEN658 IF( (nbondi == 1 .OR. nbondi == 2) .AND. (.NOT. trimmed(eidx,narea) ) ) THEN 652 659 IF( jperio /= 1 ) mbathy(nlci,:) = 0 653 660 ENDIF … … 1188 1195 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 1189 1196 USE wrk_nemo, ONLY: zenv => wrk_2d_1 , ztmp => wrk_2d_2 , zmsk => wrk_2d_3 1190 USE wrk_nemo, ONLY: z ri => wrk_2d_4 , zrj => wrk_2d_5 , zhbat => wrk_2d_61197 USE wrk_nemo, ONLY: ztmp2 => wrk_2d_4 , zhbat => wrk_2d_5 1191 1198 USE wrk_nemo, ONLY: gsigw3 => wrk_3d_1 1192 1199 USE wrk_nemo, ONLY: gsigt3 => wrk_3d_2 … … 1199 1206 USE wrk_nemo, ONLY: esigwu3 => wrk_3d_9 1200 1207 USE wrk_nemo, ONLY: esigwv3 => wrk_3d_10 1208 USE mapcomm_mod, ONLY: trimmed, cyclic_bc 1209 USE mapcomm_mod, ONLY: nidx, eidx, sidx, widx 1210 ! USE arpdebugging, ONLY: dump_array 1201 1211 !! DCSE_NEMO: wrk_nemo module variables renamed, need additional directives 1202 1212 !FTRANS gsigw3 :I :I :z … … 1213 1223 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1214 1224 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 1215 REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper ! temporary scalars 1225 REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper, zri, zrj ! temporary scalars 1226 REAL(wp), PARAMETER :: TOL_ZERO = 1.0E-20_wp ! Any value less than this assumed zero 1216 1227 ! 1217 1228 … … 1219 1230 !!---------------------------------------------------------------------- 1220 1231 1221 IF( wrk_in_use(2, 1,2,3,4,5 ,6) .OR. wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10) ) THEN1232 IF( wrk_in_use(2, 1,2,3,4,5) .OR. wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10) ) THEN 1222 1233 CALL ctl_stop('zgr_sco: ERROR - requested workspace arrays unavailable') ; RETURN 1223 1234 ENDIF … … 1269 1280 END DO 1270 1281 END DO 1282 1283 CALL lbc_lnk( zenv, 'T', 1._wp, lzero=.FALSE. ) 1271 1284 ! 1272 1285 ! Smooth the bathymetry (if required) … … 1282 1295 zrmax = 0._wp 1283 1296 zmsk(:,:) = 0._wp 1297 1284 1298 DO jj = 1, nlcj 1285 1299 DO ji = 1, nlci 1286 1300 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 1287 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 1288 zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1289 zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1290 zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 1291 IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1292 IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1._wp 1293 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1294 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1._wp 1295 END DO 1296 END DO 1301 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last row (jj=nclj+1 to jpj) 1302 zri = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1303 zrj = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1304 zrmax = MAX( zrmax, zri, zrj ) 1305 IF( zri > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1306 IF( zri > rn_rmax ) zmsk(iip1,jj ) = 1._wp 1307 IF( zrj > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1308 IF( zrj > rn_rmax ) zmsk(ji ,ijp1) = 1._wp 1309 END DO 1310 END DO 1311 1312 ! lateral boundary condition on zmsk: retain any 1's along closed 1313 ! boundary (use of lzero flag to lbc_lnk) 1314 CALL lbc_lnk( zmsk, 'T', 1._wp, lzero=.FALSE. ) 1315 1297 1316 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 1298 ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)1299 ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1._wp )1300 DO jj = 1, nlcj1301 DO ji = 1, nlci1302 zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )1303 END DO1304 END DO1305 1317 ! 1306 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 1318 IF(lwp)WRITE(numout,"('zgr_sco : iter=',I5,' rmax=',F8.4,' nb of pt= ',I8)") & 1319 jl, zrmax, INT( SUM(zmsk(:,:) ) ) 1307 1320 ! 1308 DO jj = 1, nlcj 1309 DO ji = 1, nlci 1321 !!$ IF(jl < 6)THEN ! .OR. (MOD(jl,1000) == 0) )THEN 1322 !!$ CALL dump_array(jl, 'zenv_before', zenv, withHalos=.TRUE.) 1323 !!$ CALL dump_array(jl, 'ztmp_before', ztmp, withHalos=.TRUE.) 1324 !!$ CALL dump_array(jl, 'zmsk_before', zmsk, withHalos=.TRUE.) 1325 !!$ END IF 1326 1327 ! Copy current surface before next smoothing iteration 1328 ztmp(:,:) = zenv(:,:) 1329 1330 DO jj = nldj, nlcj 1331 DO ji = nldi, nlci 1310 1332 iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci) 1311 1333 ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj) … … 1326 1348 END DO 1327 1349 ! 1328 DO jj = 1, nlcj 1329 DO ji = 1, nlci 1330 IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 1350 1351 ! Need to update halos of ztmp here but do not zero halos on closed 1352 ! boundaries 1353 CALL lbc_lnk( ztmp, 'T', 1._wp, lzero=.FALSE.) 1354 1355 DO jj = 1,nlcj 1356 DO ji = 1,nlci 1357 IF( zmsk(ji,jj) >= 1._wp-TOL_ZERO ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 1331 1358 END DO 1332 1359 END DO 1333 1360 ! 1334 ! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero 1335 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp ) 1336 DO jj = 1, nlcj 1337 DO ji = 1, nlci 1338 IF( zenv(ji,jj) == 0._wp ) zenv(ji,jj) = ztmp(ji,jj) 1339 END DO 1340 END DO 1361 ! Apply lateral boundary condition but do not zero on closed boundaries 1362 CALL lbc_lnk( zenv, 'T', 1._wp, lzero=.FALSE. ) 1363 1364 !!$ IF(jl < 6)THEN ! .OR. (MOD(jl,1000) == 0) )THEN 1365 !!$ CALL dump_array(jl, 'zenv', zenv, withHalos=.TRUE.) 1366 !!$ CALL dump_array(jl, 'ztmp', ztmp, withHalos=.TRUE.) 1367 !!$ CALL dump_array(jl, 'zmsk', zmsk, withHalos=.TRUE.) 1368 !!$ END IF 1369 1341 1370 ! ! ================ ! 1342 1371 END DO ! End loop ! … … 1365 1394 ENDIF 1366 1395 ENDIF 1396 1397 ! CALL dump_array(0, 'hbatt', hbatt, withHalos=.FALSE.) 1367 1398 1368 1399 ! ! ============================== … … 1576 1607 ! 1577 1608 !! H. Liu, POL. April 2009. Added for passing the scale check for the new released vvl code. 1609 where (e3t (:,:,:).eq.0.0) e3t(:,:,:) = 1.0 1610 where (e3u (:,:,:).eq.0.0) e3u(:,:,:) = 1.0 1611 where (e3v (:,:,:).eq.0.0) e3v(:,:,:) = 1.0 1612 where (e3f (:,:,:).eq.0.0) e3f(:,:,:) = 1.0 1613 where (e3w (:,:,:).eq.0.0) e3w(:,:,:) = 1.0 1614 where (e3uw (:,:,:).eq.0.0) e3uw(:,:,:) = 1.0 1615 where (e3vw (:,:,:).eq.0.0) e3vw(:,:,:) = 1.0 1616 1578 1617 1579 1618 fsdept(:,:,:) = gdept (:,:,:) … … 1597 1636 END DO 1598 1637 END DO 1638 1599 1639 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & 1600 1640 & ' MAX ', MAXVAL( mbathy(:,:) ) … … 1642 1682 END DO 1643 1683 END DO 1644 DO jj = mj0(74), mj1(74) 1645 DO ji = mi0(100), mi1(100) 1646 WRITE(numout,*) 1647 WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1648 WRITE(numout,*) ' ~~~~~~ --------------------' 1649 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1650 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & 1651 & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 1652 END DO 1653 END DO 1684 !!$ ARPDBG - out of bounds if jpj < 74 or jpi < 100, e.g. default GYRE 1685 !!$ DO jj = mj0(74), mj1(74) 1686 !!$ DO ji = mi0(100), mi1(100) 1687 !!$ WRITE(numout,*) 1688 !!$ WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1689 !!$ WRITE(numout,*) ' ~~~~~~ --------------------' 1690 !!$ WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1691 !!$ WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & 1692 !!$ & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 1693 !!$ END DO 1694 !!$ END DO 1654 1695 ENDIF 1655 1696 … … 1668 1709 CALL ctl_stop( ctmp1 ) 1669 1710 ENDIF 1670 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN1711 IF( gdepw_1(ji,jj,jk) < 0._wp .OR. gdept_1(ji,jj,jk) < 0._wp ) THEN 1671 1712 WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1672 1713 CALL ctl_stop( ctmp1 ) … … 1677 1718 !!gm bug #endif 1678 1719 ! 1679 IF( wrk_not_released(2, 1,2,3,4,5 ,6) .OR. wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10) ) &1720 IF( wrk_not_released(2, 1,2,3,4,5) .OR. wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10) ) & 1680 1721 & CALL ctl_stop('dom:zgr_sco: failed to release workspace arrays') 1681 1722 ! -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r3211 r3432 144 144 IF( lk_vvl ) THEN 145 145 #if defined key_z_first 146 fse3t_b(:,:,:) = fse3t_n(:,:,:) 146 ! DCSE_NEMO: can't use implicit loop over k here because the domzgr_substitute.h90 147 ! file causes the line below to be expanded to: 148 ! e3t_b(:,:,:) = (e3t(:,:,:)*(1.+sshn(:,:)*mut(:,:,:))) 149 ! which contains non-conforming array expressions. 150 DO jk = 1, jpk 151 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 152 END DO 147 153 #else 148 154 DO jk = 1, jpk … … 501 507 502 508 IF(wrk_in_use(3, 1) ) THEN 503 CALL ctl_stop('ista ge_uvg: requested workspace array unavailable') ; RETURN509 CALL ctl_stop('istate_uvg: requested workspace array unavailable') ; RETURN 504 510 ENDIF 505 511 … … 628 634 ! 629 635 IF( wrk_not_released(3, 1) ) THEN 630 CALL ctl_stop('ista ge_uvg: failed to release workspace array')636 CALL ctl_stop('istate_uvg: failed to release workspace array') 631 637 ENDIF 632 638 ! -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r3211 r3432 70 70 71 71 IF( wrk_in_use(3,1) ) THEN 72 CALL ctl_stop('dyn_ke y: requested workspace array is unavailable') ; RETURN72 CALL ctl_stop('dyn_keg: requested workspace array is unavailable') ; RETURN 73 73 ENDIF 74 74 … … 164 164 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 165 165 ! 166 IF( wrk_not_released(3, 1) ) CALL ctl_stop('dyn_ke y: failed to release workspace array')166 IF( wrk_not_released(3, 1) ) CALL ctl_stop('dyn_keg: failed to release workspace array') 167 167 ! 168 168 END SUBROUTINE dyn_keg -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
r3211 r3432 33 33 PUBLIC dyn_ldf_bilapg ! called by step.F90 34 34 35 !FTRANS zfuw zfvw zdiu zdiv :I :z 35 36 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zfvw , zdiu, zdiv ! 2D workspace (ldfguv) 37 !FTRANS zdju zdj1u zdjv zdj1v :I :z 36 38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdju, zdj1u, zdjv, zdj1v ! 2D workspace (ldfguv) 37 39 … … 113 115 IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate' 114 116 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 115 zwk1(:,:,:) = 0.e0 ; zwk3(:,:,:) = 0.e0116 zwk2(:,:,:) = 0.e0 ; zwk4(:,:,:) = 0.e0117 zwk1(:,:,:) = 0.e0_wp ; zwk3(:,:,:) = 0.e0_wp 118 zwk2(:,:,:) = 0.e0_wp ; zwk4(:,:,:) = 0.e0_wp 117 119 ! ! allocate dyn_ldf_bilapg arrays 118 120 IF( dyn_ldf_bilapg_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays') … … 195 197 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 196 198 USE wrk_nemo, ONLY: ziut => wrk_2d_1 , zjuf => wrk_2d_2 , zjvt => wrk_2d_3 197 USE wrk_nemo, ONLY: zivf => wrk_2d_4 , zdku => wrk_2d_5 , zdk1u => wrk_2d_6 199 USE wrk_nemo, ONLY: zivf => wrk_2d_4 200 #if ! defined key_z_first 201 USE wrk_nemo, ONLY: zdku => wrk_2d_5 , zdk1u => wrk_2d_6 198 202 USE wrk_nemo, ONLY: zdkv => wrk_2d_7 , zdk1v => wrk_2d_8 203 #endif 204 USE timing , ONLY: timing_start, timing_stop 199 205 !! 200 206 !FTRANS pu :I :I :z … … 213 219 ! 214 220 INTEGER :: ji, jj, jk ! dummy loop indices 221 INTEGER :: jif, jjf ! dummy loop indices over full domain 215 222 REAL(wp) :: zabe1 , zabe2 , zcof1 , zcof2 ! local scalar 216 223 REAL(wp) :: zcoef0, zcoef3, zcoef4 ! - - 217 224 REAL(wp) :: zbur, zbvr, zmkt, zmkf, zuav, zvav ! - - 218 225 REAL(wp) :: zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - - 226 #if defined key_z_first 227 ! Can use scalars instead of work arrays when built with z-first 228 REAL(wp) :: zdku, zdkv, zdk1u, zdk1v 229 #endif 219 230 !!---------------------------------------------------------------------- 231 232 CALL timing_start('ldfguv') 220 233 221 234 IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN 222 235 CALL ctl_stop('dyn:ldfguv: requested workspace arrays unavailable') ; RETURN 223 236 END IF 237 CALL timing_start('ldfguv_1st') 238 239 #if defined key_z_first 240 ! ! ********** ! ! =============== 241 ! DO jk = 1, jpkm1 ! First step ! ! Horizontal slab 242 ! ! ********** ! ! =============== 243 244 ! I.1 Vertical gradient of pu and pv at level jk and jk+1 245 ! ------------------------------------------------------- 246 ! surface boundary condition: zdku(jk=1)=zdku(jk=2) 247 ! zdkv(jk=1)=zdkv(jk=2) 248 #if 0 249 DO jjf = 1, jpj 250 DO jif = 1, jpi 251 252 !!$ jj= jjf 253 !!$ ji = jif 254 !!$ 255 !!$ zdk1u = ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) 256 !!$ zdk1v = ( pv(ji,jj,jk) - pv(ji,jj,jk+1) ) * vmask(ji,jj,jk+1) 257 !!$ 258 !!$ IF( jk == 1 ) THEN 259 !!$ zdku = zdk1u 260 !!$ zdkv = zdk1v 261 !!$ ELSE 262 !!$ zdku = ( pu(ji,jj,jk-1) - pu(ji,jj,jk) ) * umask(ji,jj,jk) 263 !!$ zdkv = ( pv(ji,jj,jk-1) - pv(ji,jj,jk) ) * vmask(ji,jj,jk) 264 !!$ ENDIF 265 ! END DO 266 ! END DO 267 ! -----f----- 268 ! I.2 Horizontal fluxes on U | 269 ! ------------------------=== t u t 270 ! | 271 ! i-flux at t-point -----f----- 272 ! DO jj = 1, jpjm1 273 ! DO ji = 2, jpi 274 ! DO jjf = 1, jpj 275 ! DO jif = 1, jpi 276 277 jj = MIN(jjf, jpjm1) 278 ji = MAX(2, jif) 279 280 ! zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 281 282 ! zmkt = 1._wp/MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 283 ! + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1._wp ) 284 285 ! zcof1 = -e2t(ji,jj) / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 286 ! + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1._wp ) & 287 ! * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 288 289 IF( jk == 1 )THEN 290 ziut(ji,jj) = tmask(ji,jj,jk) * & 291 ( (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & 292 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 293 + (-e2t(ji,jj) / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 294 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1._wp ) & 295 * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * & 296 ( ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + & 297 ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + & 298 ! +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) 299 ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + & 300 ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) & 301 ) ) 302 ELSE 303 ziut(ji,jj) = tmask(ji,jj,jk) * & 304 ( (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & 305 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 306 + (-e2t(ji,jj) / MAX(umask(ji-1,jj,jk)+umask(ji,jj,jk+1) & 307 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk), 1._wp ) & 308 * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * & 309 ( ( pu(ji ,jj,jk-1) - pu(ji ,jj,jk ) ) * umask(ji,jj,jk) + & 310 ( pu(ji-1,jj,jk ) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + & 311 ! +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) 312 ( pu(ji ,jj,jk ) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + & 313 ( pu(ji-1,jj,jk-1) - pu(ji-1,jj,jk ) ) * umask(ji-1,jj,jk) & 314 ) ) 315 END IF 316 ! END DO 317 ! END DO 318 319 ! j-flux at f-point 320 ! DO jj = 1, jpjm1 321 ! DO ji = 1, jpim1 322 ! DO jjf = 1, jpj 323 ! DO jif = 1, jpi 324 325 jj = MIN(jjf, jpjm1) 326 ji = MIN(jif, jpim1) 327 328 !zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 329 330 !zmkf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & 331 ! + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. ) 332 333 !zcof2 = -e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1) & 334 ! + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. ) & 335 ! * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 336 337 IF(jk == 1)THEN 338 ! zdku = zdk1u 339 zjuf(ji,jj) = fmask(ji,jj,jk) * & 340 ( (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & 341 ! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & 342 + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1) & 343 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. ) & 344 * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * & 345 ( & 346 ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) & 347 + ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj ,jk+1) & 348 ! +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) 349 + ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) & 350 + ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) & 351 ) ) 352 353 ELSE 354 zjuf(ji,jj) = fmask(ji,jj,jk) * & 355 ( (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & 356 ! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & 357 + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1) & 358 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. ) & 359 * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * & 360 ( & 361 (pu(ji,jj+1,jk-1) - pu(ji,jj+1,jk ) ) * umask(ji,jj+1,jk) + & 362 (pu(ji,jj ,jk ) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + & 363 ! +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) 364 (pu(ji,jj+1,jk ) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) + & 365 (pu(ji,jj ,jk-1) - pu(ji,jj ,jk ) ) * umask(ji,jj,jk) & 366 ) ) 367 ENDIF 368 ! END DO 369 ! END DO 370 371 ! | t | 372 ! I.3 Horizontal fluxes on V | | 373 ! ------------------------=== f---v---f 374 ! | | 375 ! i-flux at f-point | t | 376 ! DO jj = 1, jpjm1 377 ! DO ji = 1, jpim1 378 ! DO jjf = 1, jpj 379 ! DO jif = 1, jpi 380 381 jj = MIN(jjf, jpjm1) 382 ji = MIN(jif, jpim1) 383 384 ! zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 385 386 ! zmkf = 1./MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 387 ! + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. ) 388 389 ! zcof1 = (-e2f(ji,jj) / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 390 ! + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. )) & 391 ! * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 392 393 IF (jk == 1)THEN 394 ! zdku == zdk1u 395 zivf(ji,jj) = fmask(ji,jj,jk) * & 396 ( (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & 397 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & 398 + ((-e2f(ji,jj) / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 399 + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. )) & 400 * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( & 401 ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & 402 ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & 403 ! +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) 404 ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & 405 ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) ) ) 406 ELSE 407 zivf(ji,jj) = fmask(ji,jj,jk) * & 408 ( (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & 409 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & 410 + ((-e2f(ji,jj) / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 411 + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. )) & 412 * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( & 413 ( pu(ji ,jj,jk-1) - pu(ji ,jj,jk ) ) * umask(ji ,jj,jk ) + & 414 ( pu(ji+1,jj,jk ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & 415 ! +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) 416 ( pu(ji ,jj,jk ) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & 417 ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk ) ) * umask(ji+1,jj,jk ) ) ) 418 END IF 419 ! END DO 420 ! END DO 421 422 ! j-flux at t-point 423 ! DO jj = 2, jpj 424 ! DO ji = 1, jpim1 425 ! DO jjf = 1, jpj 426 ! DO jif = 1, jpi 427 428 jj = MAX(2,jjf) 429 ji = MIN(jif, jpim1) 430 431 !zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 432 433 !zmkt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 434 ! + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) 435 436 !zcof2 = (-e1t(ji,jj)/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 437 ! + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) ) & 438 ! * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 439 440 IF( jk == 1 )THEN 441 zjvt(ji,jj) = tmask(ji,jj,jk) * & 442 ( (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & 443 ! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & 444 + ((-e1t(ji,jj)/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 445 + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) ) & 446 * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * ( & 447 ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + & 448 ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + & 449 ! +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) 450 ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + & 451 ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) ) ) 452 ELSE 453 zjvt(ji,jj) = tmask(ji,jj,jk) * & 454 ( (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & 455 ! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & 456 + ((-e1t(ji,jj)/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 457 + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) ) & 458 * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * & 459 ( & 460 ( pu(ji,jj-1,jk-1) - pu(ji,jj-1,jk ) ) * umask(ji,jj-1,jk) + & 461 ( pu(ji,jj ,jk ) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + & 462 ! +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) 463 ( pu(ji,jj-1,jk ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + & 464 ( pu(ji,jj ,jk-1) - pu(ji,jj ,jk ) ) * umask(ji,jj,jk) ) ) 465 END IF 466 END DO 467 END DO 468 #endif 469 470 ! I.4 Second derivative (divergence) (not divided by the volume) 471 ! --------------------- 472 473 474 DO jj = 2, jpjm1 475 DO ji = 2, jpim1 476 477 ! Treat jk = 1 separately as is special case 478 jk = 1 479 480 plu(ji,jj,jk) = & 481 ! ------------- ziut (ji+1, jj) - 482 tmask(ji+1,jj,jk) * & 483 ( (e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & 484 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 485 + (-e2t(ji+1,jj) / MAX( umask(ji,jj,jk )+umask(ji+1,jj,jk+1) & 486 + umask(ji,jj,jk+1)+umask(ji+1,jj,jk ), 1._wp ) & 487 * 0.5 * ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) ) ) * & 488 ( ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & 489 ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & 490 ! +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) 491 ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & 492 ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) & 493 ) ) - & 494 ! ------------- ziut (ji,jj ) + 495 tmask(ji,jj,jk) * & 496 ( (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & 497 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 498 + (-e2t(ji,jj) / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 499 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1._wp ) & 500 * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * & 501 ( ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + & 502 ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + & 503 ! +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) 504 ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + & 505 ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) & 506 ) ) + & 507 ! ------------- zjuf (ji ,jj) - 508 fmask(ji,jj,jk) * & 509 ( (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & 510 ! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & 511 + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1) & 512 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. ) & 513 * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * & 514 ( & 515 ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) & 516 + ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj ,jk+1) & 517 ! +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) 518 + ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) & 519 + ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) & 520 ) ) - & 521 ! ------------- zjuf (ji,jj-1) 522 fmask(ji,jj-1,jk) * & 523 ( (e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1) ) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & 524 ! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & 525 + (-e1f(ji,jj-1) /MAX(umask(ji,jj,jk)+umask(ji,jj-1,jk+1) & 526 + umask(ji,jj,jk+1)+umask(ji,jj-1,jk), 1. ) & 527 * 0.5 * ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) ) ) * & 528 ( & 529 ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) & 530 + ( pu(ji,jj-1 ,jk) - pu(ji,jj-1 ,jk+1) ) * umask(ji,jj-1 ,jk+1) & 531 ! +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) 532 + ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) & 533 + ( pu(ji,jj-1 ,jk) - pu(ji,jj-1 ,jk+1) ) * umask(ji,jj-1,jk+1) & 534 ) ) 535 536 537 plv(ji,jj,jk) = & 538 ! ------------- zivf (ji,jj ) - 539 fmask(ji,jj,jk) * & 540 ( (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & 541 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & 542 + ((-e2f(ji,jj) / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 543 + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. )) & 544 * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( & 545 ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & 546 ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & 547 ! +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) 548 ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & 549 ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) & 550 ) ) - & 551 ! ------------- zivf (ji-1,jj) + 552 fmask(ji-1,jj,jk) * & 553 ( (e2f(ji-1,jj) * fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & 554 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & 555 + ((-e2f(ji-1,jj) / MAX( vmask(ji,jj,jk )+vmask(ji-1,jj,jk+1) & 556 + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk ), 1. )) & 557 * 0.5 * ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * ( & 558 ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + & 559 ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & 560 ! +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) 561 ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + & 562 ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) & 563 ) ) + & 564 ! ------------- zjvt (ji,jj+1) - 565 tmask(ji,jj+1,jk) * & 566 ( (e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / e2t(ji,jj+1)) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & 567 ! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & 568 + ((-e1t(ji,jj+1)/MAX( vmask(ji,jj,jk )+vmask(ji,jj+1,jk+1) & 569 + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk ), 1. ) ) & 570 * 0.5 * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * ( & 571 ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj ,jk+1) + & 572 ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) + & 573 ! +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) 574 ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj ,jk+1) + & 575 ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) & 576 ) ) - & 577 ! ------------- zjvt (ji,jj ) 578 tmask(ji,jj,jk) * & 579 ( (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & 580 ! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & 581 + ((-e1t(ji,jj)/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 582 + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) ) & 583 * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * ( & 584 ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + & 585 ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + & 586 ! +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) 587 ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + & 588 ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) & 589 ) ) 590 591 592 DO jk = 2, jpkm1 593 594 ! plu(ji,jj,jk) = ziut (ji+1,jj) - & 595 ! ziut (ji,jj ) + & 596 ! zjuf (ji ,jj) - & 597 ! zjuf (ji,jj-1) 598 plu(ji,jj,jk) = & 599 ! ------------- ziut (ji+1, jj ) - 600 tmask(ji+1,jj,jk) * & 601 ( (e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & 602 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 603 + (-e2t(ji+1,jj) / MAX(umask(ji,jj,jk)+umask(ji+1,jj,jk+1) & 604 + umask(ji,jj,jk+1)+umask(ji+1,jj,jk), 1._wp ) & 605 * 0.5 * ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) ) ) * & 606 ( ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk ) ) * umask(ji+1,jj,jk) + & 607 ( pu(ji ,jj,jk ) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & 608 ! +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) 609 ( pu(ji+1,jj,jk ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & 610 ( pu(ji ,jj,jk-1) - pu(ji ,jj,jk ) ) * umask(ji ,jj,jk) & 611 ) ) - & 612 ! ------------- ziut (ji , jj ) + 613 tmask(ji,jj,jk) * & 614 ( (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & 615 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 616 + (-e2t(ji,jj) / MAX(umask(ji-1,jj,jk)+umask(ji,jj,jk+1) & 617 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk), 1._wp ) & 618 * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * & 619 ( ( pu(ji ,jj,jk-1) - pu(ji ,jj,jk ) ) * umask(ji,jj,jk) + & 620 ( pu(ji-1,jj,jk ) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + & 621 ! +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) 622 ( pu(ji ,jj,jk ) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + & 623 ( pu(ji-1,jj,jk-1) - pu(ji-1,jj,jk ) ) * umask(ji-1,jj,jk) & 624 ) ) + & 625 ! ------------- zjuf (ji , jj ) - 626 fmask(ji,jj,jk) * & 627 ( (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & 628 ! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & 629 + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1) & 630 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. ) & 631 * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * & 632 ( & 633 (pu(ji,jj+1,jk-1) - pu(ji,jj+1,jk ) ) * umask(ji,jj+1,jk) + & 634 (pu(ji,jj ,jk ) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + & 635 ! +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) 636 (pu(ji,jj+1,jk ) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) + & 637 (pu(ji,jj ,jk-1) - pu(ji,jj ,jk ) ) * umask(ji,jj,jk) & 638 ) ) - & 639 ! ------------- zjuf (ji , jj-1) 640 fmask(ji,jj-1,jk) * & 641 ( (e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1) ) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & 642 ! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & 643 + (-e1f(ji,jj-1) /MAX(umask(ji,jj,jk)+umask(ji,jj-1,jk+1) & 644 + umask(ji,jj,jk+1)+umask(ji,jj-1,jk), 1. ) & 645 * 0.5 * ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) ) ) * & 646 ( & 647 (pu(ji,jj,jk-1) - pu(ji,jj,jk ) ) * umask(ji,jj,jk) + & 648 (pu(ji,jj-1,jk ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + & 649 ! +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) 650 (pu(ji,jj,jk ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + & 651 (pu(ji,jj-1,jk-1) - pu(ji,jj-1 ,jk ) ) * umask(ji,jj-1,jk) & 652 ) ) 653 654 655 ! plv(ji,jj,jk) = zivf (ji,jj ) - & 656 ! zivf (ji-1,jj) + & 657 ! zjvt (ji,jj+1) - & 658 ! zjvt (ji,jj ) 659 plv(ji,jj,jk) = & 660 ! ------------- zivf (ji,jj ) - 661 fmask(ji,jj,jk) * & 662 ( (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & 663 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & 664 + ((-e2f(ji,jj) / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 665 + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. )) & 666 * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( & 667 ( pu(ji ,jj,jk-1) - pu(ji ,jj,jk ) ) * umask(ji ,jj,jk ) + & 668 ( pu(ji+1,jj,jk ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & 669 ! +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) 670 ( pu(ji ,jj,jk ) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & 671 ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk ) ) * umask(ji+1,jj,jk ) & 672 ) ) - & 673 ! ------------- zivf (ji-1,jj) + 674 fmask(ji-1,jj,jk) * & 675 ( (e2f(ji-1,jj) * fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & 676 ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & 677 + ((-e2f(ji-1,jj) / MAX( vmask(ji,jj,jk )+vmask(ji-1,jj,jk+1) & 678 + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk ), 1. )) & 679 * 0.5 * ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * ( & 680 ( pu(ji-1 ,jj,jk-1) - pu(ji-1 ,jj,jk ) ) * umask(ji-1 ,jj,jk ) + & 681 ( pu(ji,jj,jk ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + & 682 ! +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) 683 ( pu(ji-1 ,jj,jk ) - pu(ji-1 ,jj,jk+1) ) * umask(ji-1 ,jj,jk+1) + & 684 ( pu(ji,jj,jk-1) - pu(ji,jj,jk ) ) * umask(ji,jj,jk ) & 685 ) ) + & 686 ! ------------- zjvt (ji,jj+1) - 687 tmask(ji,jj+1,jk) * & 688 ( (e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / e2t(ji,jj+1)) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & 689 ! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & 690 + ((-e1t(ji,jj+1)/MAX( vmask(ji,jj,jk )+vmask(ji,jj+1,jk+1) & 691 + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk ), 1. ) ) & 692 * 0.5 * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * & 693 ( & 694 ( pu(ji,jj,jk-1) - pu(ji,jj,jk ) ) * umask(ji,jj,jk) + & 695 ( pu(ji,jj+1 ,jk ) - pu(ji,jj+1 ,jk+1) ) * umask(ji,jj+1,jk+1) + & 696 ! +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) 697 ( pu(ji,jj,jk ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + & 698 ( pu(ji,jj+1 ,jk-1) - pu(ji,jj+1 ,jk ) ) * umask(ji,jj+1,jk) & 699 ) ) - & 700 ! ------------- zjvt (ji,jj ) 701 tmask(ji,jj,jk) * & 702 ( (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & 703 ! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & 704 + ((-e1t(ji,jj)/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 705 + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) ) & 706 * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * & 707 ( & 708 ( pu(ji,jj-1,jk-1) - pu(ji,jj-1,jk ) ) * umask(ji,jj-1,jk) + & 709 ( pu(ji,jj ,jk ) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + & 710 ! +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) 711 ( pu(ji,jj-1,jk ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + & 712 ( pu(ji,jj ,jk-1) - pu(ji,jj ,jk ) ) * umask(ji,jj,jk) ) ) 713 714 END DO 715 END DO 716 717 ! ! =============== 718 END DO ! End of slab 719 ! ! =============== 720 #else 224 721 ! ! ********** ! ! =============== 225 722 DO jk = 1, jpkm1 ! First step ! ! Horizontal slab … … 338 835 END DO ! End of slab 339 836 ! ! =============== 837 #endif 838 CALL timing_stop('ldfguv_1st','section') 340 839 341 840 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 342 841 842 CALL timing_start('ldfguv_2nd') 343 843 ! ! ************ ! ! =============== 344 844 DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab … … 348 848 ! --------------------------------- 349 849 850 #if defined key_z_first 851 DO ji = 2, jpi 852 DO jk = 1, jpk 853 ! i-gradient of u at jj 854 zdiu (ji,jk) = tmask(ji,jj ,jk) * ( pu(ji,jj ,jk) - pu(ji-1,jj ,jk) ) 855 ! j-gradient of u and v at jj 856 zdju (ji,jk) = fmask(ji,jj ,jk) * ( pu(ji,jj+1,jk) - pu(ji ,jj ,jk) ) 857 zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pv(ji,jj ,jk) - pv(ji ,jj-1,jk) ) 858 ! j-gradient of u and v at jj+1 859 zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj ,jk) - pu(ji ,jj-1,jk) ) 860 zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji ,jj ,jk) ) 861 END DO 862 END DO 863 DO ji = 1, jpim1 864 DO jk = 1, jpk 865 ! i-gradient of v at jj 866 zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pv(ji+1,jj,jk) - pv(ji ,jj ,jk) ) 867 END DO 868 END DO 869 #else 350 870 DO jk = 1, jpk 351 871 DO ji = 2, jpi … … 366 886 END DO 367 887 END DO 368 888 #endif 369 889 370 890 ! II.2 Vertical fluxes … … 380 900 ! interior (2=<jk=<jpk-1) on pu field 381 901 902 #if defined key_z_first 903 DO ji = 2, jpim1 904 DO jk = 2, jpkm1 905 #else 382 906 DO jk = 2, jpkm1 383 907 DO ji = 2, jpim1 908 #endif 384 909 ! i- and j-slopes at uw-point 385 910 zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) … … 408 933 ! interior (2=<jk=<jpk-1) on pv field 409 934 935 #if defined key_z_first 936 DO ji = 2, jpim1 937 DO jk = 2, jpkm1 938 #else 410 939 DO jk = 2, jpkm1 411 940 DO ji = 2, jpim1 941 #endif 412 942 ! i- and j-slopes at vw-point 413 943 zvwslpi = 0.5 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) … … 440 970 IF( kahm == 1 ) THEN 441 971 ! multiply the laplacian by the eddy viscosity coefficient 972 #if defined key_z_first 973 DO ji = 2, jpim1 974 DO jk = 1, jpkm1 975 #else 442 976 DO jk = 1, jpkm1 443 977 DO ji = 2, jpim1 978 #endif 444 979 ! eddy coef. divided by the volume element 445 980 zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) ) … … 455 990 ELSEIF( kahm == 2 ) THEN 456 991 ! second call, no multiplication 992 #if defined key_z_first 993 DO ji = 2, jpim1 994 DO jk = 1, jpkm1 995 #else 457 996 DO jk = 1, jpkm1 458 997 DO ji = 2, jpim1 998 #endif 459 999 ! inverse of the volume element 460 1000 zbur = 1. / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) ) … … 476 1016 END DO ! End of slab 477 1017 ! ! =============== 1018 CALL timing_stop('ldfguv_2nd','section') 478 1019 479 1020 IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) ) CALL ctl_stop('dyn:ldfguv: failed to release workspace arrays') 480 1021 ! 1022 CALL timing_stop('ldfguv','section') 1023 481 1024 END SUBROUTINE ldfguv 482 1025 -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r3211 r3432 36 36 PUBLIC dyn_ldf_iso_alloc ! called by nemogcm.F90 37 37 38 !FTRANS zdiu zdju zdiv zdjv zdj1u zdj1v zfuw zfvw :I :z 38 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso) 39 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v ! - - … … 114 115 !!---------------------------------------------------------------------- 115 116 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 116 USE wrk_nemo, ONLY: ziut => wrk_2d_1 , zjuf => wrk_2d_2 , zjvt => wrk_2d_3 ! 2D workspace 117 USE wrk_nemo, ONLY: zivf => wrk_2d_4 , zdku => wrk_2d_5 , zdkv => wrk_2d_6 ! 2D workspace 117 #if ! defined key_z_first 118 ! Then these workspace arrays can be left as 2D 119 USE wrk_nemo, ONLY: zjvt => wrk_2d_3 ! 2D workspace 120 USE wrk_nemo, ONLY: zivf => wrk_2d_4 ! 2D workspace 121 USE wrk_nemo, ONLY: ziut => wrk_2d_1 , zjuf => wrk_2d_2 122 USE wrk_nemo, ONLY: zdku => wrk_2d_5 , zdkv => wrk_2d_6 118 123 USE wrk_nemo, ONLY: zdk1u => wrk_2d_7 , zdk1v => wrk_2d_8 124 #endif 125 USE timing, ONLY: timing_start, timing_stop 119 126 ! 120 127 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 125 132 REAL(wp) :: zcoef0, zcoef3, zcoef4, zmkt, zmkf ! - - 126 133 REAL(wp) :: zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - - 134 REAL(wp) :: zcof, zrecip 135 #if defined key_z_first 136 REAL(wp) :: zdku, zdk1u, zdki1u, zdk1i1u, zdkim1u, zdk1im1u 137 REAL(wp) :: zdkj1u, zdk1j1u, zdkjm1u, zdk1jm1u 138 REAL(wp) :: zdkv, zdk1v, zdki1v, zdk1i1v, zdkim1v, zdk1im1v 139 REAL(wp) :: zdkj1v, zdk1j1v, zdkjm1v, zdk1jm1v 140 REAL(wp) :: aziut, aziuti1, azjuf, azjufjm1 141 REAL(wp) :: azivf, azivfim1, azjvtj1, azjvt 142 #endif 127 143 !!---------------------------------------------------------------------- 144 145 ! ziut zjvt zjuf zivf :I :I :z 146 !!$#if defined key_z_first 147 !!$ REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ziut, zjuf, zivf, zjvt 148 !!$ ALLOCATE( ziut(jpi,jpj,jpk), zjuf(jpi,jpj,jpk), & 149 !!$ zivf(jpi,jpj,jpk), zjvt(jpi,jpj,jpk) ) 150 !!$#endif 151 CALL timing_start('dyn_ldf_iso') 128 152 129 153 IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN … … 170 194 ENDIF 171 195 196 CALL timing_start('dyn_ldf_iso_hslab') 197 198 #if defined key_z_first 199 200 ! Vertical u- and v-shears at level jk and jk+1 201 ! --------------------------------------------- 202 ! surface boundary condition: zdku(jk=1)=zdku(jk=2) 203 ! zdkv(jk=1)=zdkv(jk=2) 204 !!$ DO jj = 1, jpj, 1 205 !!$ DO ji = 1, jpi, 1 206 !!$ 207 !!$ ! jk=1 special case 208 !!$ !zdk1u(ji,jj,1) = ( ub(ji,jj,1) -ub(ji,jj,2) ) * umask(ji,jj,2) 209 !!$ zdk1v(ji,jj,1) = ( vb(ji,jj,1) -vb(ji,jj,2) ) * vmask(ji,jj,2) 210 !!$ !zdku(ji,jj,1) = zdk1u(ji,jj,1) 211 !!$ zdkv(ji,jj,1) = zdk1v(ji,jj,1) 212 !!$ 213 !!$ DO jk = 2, jpkm1, 1 214 !!$ !zdk1u(ji,jj,jk) = ( ub(ji,jj,jk) -ub(ji,jj,jk+1) ) * umask(ji,jj,jk+1) 215 !!$ zdk1v(ji,jj,jk) = ( vb(ji,jj,jk) -vb(ji,jj,jk+1) ) * vmask(ji,jj,jk+1) 216 !!$ 217 !!$ !zdku(ji,jj,jk) = ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) * umask(ji,jj,jk) 218 !!$ zdkv(ji,jj,jk) = ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) * vmask(ji,jj,jk) 219 !!$ END DO !jk 220 !!$ END DO 221 !!$ END DO 222 223 !!$ ! -----f----- 224 !!$ ! Horizontal fluxes on U | 225 !!$ ! --------------------=== t u t 226 !!$ ! | 227 !!$ ! i-flux at t-point -----f----- 228 !!$ 229 !!$ IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) 230 !!$ DO jj = 2, jpjm1 231 !!$ DO ji = 2, jpi 232 !!$ DO jk = 1, jpkm1, 1 233 !!$ zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj) 234 !!$ 235 !!$ zmskt = 1._wp/MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 236 !!$ & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1._wp ) 237 !!$ 238 !!$ zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5_wp * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 239 !!$ 240 !!$ ziut(ji,jj,jk) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & 241 !!$ & + zcof1 * ( zdku + zdk1im1u & 242 !!$ & +zdk1u + zdkim1u ) ) * tmask(ji,jj,jk) 243 !!$ END DO 244 !!$ END DO 245 !!$ END DO 246 !!$ ELSE ! other coordinate system (zco or sco) : e3t 247 !!$ DO jj = 2, jpjm1 248 !!$ DO ji = 2, jpi 249 !!$ DO jk = 1, jpkm1, 1 250 !!$ zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 251 !!$ 252 !!$ zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 253 !!$ & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. ) 254 !!$ 255 !!$ zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 256 !!$ 257 !!$ ziut(ji,jj,jk) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & 258 !!$ & + zcof1 * ( zdku(ji,jj,jk) + zdk1u(ji-1,jj,jk) & 259 !!$ & +zdk1u(ji,jj,jk) + zdku (ji-1,jj,jk) ) ) * tmask(ji,jj,jk) 260 !!$ 261 !!$ END DO 262 !!$ END DO 263 !!$ END DO 264 !!$ ENDIF 265 266 ! j-flux at f-point 267 ! BLOCKABLE(ji,jj,jk) 268 ! BLOCKING SIZE (0) 269 !!$ DO jj = 1, jpjm1, 1 270 !!$! BLOCKING SIZE (0) 271 !!$ DO ji = 1, jpim1, 1 272 !!$! BLOCKING SIZE (4) 273 !!$ DO jk = 1, jpkm1, 1 274 !!$ zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 275 !!$ 276 !!$ zmskf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & 277 !!$ & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. ) 278 !!$ 279 !!$ zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 280 !!$ 281 !!$ zjuf(ji,jj,jk) = ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) & 282 !!$ & + zcof2 * ( zdku (ji,jj+1,jk) + zdk1u(ji,jj,jk) & 283 !!$ & +zdk1u(ji,jj+1,jk) + zdku (ji,jj,jk) ) ) * fmask(ji,jj,jk) 284 !!$ END DO 285 !!$ END DO 286 !!$ END DO 287 288 ! | t | 289 ! Horizontal fluxes on V | | 290 ! --------------------=== f---v---f 291 ! | | 292 ! | t | 293 294 !!$ IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) 295 !!$ DO jj = 2, jpj 296 !!$ DO ji = 1, jpim1 297 !!$ DO jk = 1, jpkm1, 1 298 !!$ 299 !!$ ! i-flux at f-point | t | 300 !!$ zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 301 !!$ 302 !!$ zmskf = 1._wp/MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 303 !!$ + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp ) 304 !!$ 305 !!$ zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 306 !!$ 307 !!$ zivf(ji,jj,jk) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & 308 !!$ + zcof1 * ( zdkv(ji,jj,jk) + zdk1v(ji+1,jj,jk) & 309 !!$ +zdk1v(ji,jj,jk) + zdkv(ji+1,jj,jk) ) ) * fmask(ji,jj,jk) 310 !!$ 311 !!$ 312 !!$ ! j-flux at t-point 313 !!$ zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj) 314 !!$ 315 !!$ zmskt = 1._wp/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 316 !!$ + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp ) 317 !!$ 318 !!$ zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5_wp * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 319 !!$ 320 !!$ zjvt(ji,jj,jk) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & 321 !!$ + zcof2 * ( zdkv(ji,jj-1,jk) + zdk1v(ji,jj,jk) & 322 !!$ +zdk1v(ji,jj-1,jk) + zdkv (ji,jj,jk) ) ) * tmask(ji,jj,jk) 323 !!$ END DO 324 !!$ END DO 325 !!$ END DO 326 !!$ ELSE ! other coordinate system (zco or sco) : e3t 327 !!$ DO jj = 2, jpj 328 !!$ DO ji = 1, jpim1 329 !!$ DO jk = 1, jpkm1 330 !!$ 331 !!$ ! i-flux at f-point 332 !!$ zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 333 !!$ 334 !!$ zmskf = 1._wp/MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 335 !!$ + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp ) 336 !!$ 337 !!$ zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 338 !!$ 339 !!$ zivf(ji,jj,jk) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & 340 !!$ + zcof1 * ( zdkv(ji,jj,jk) + zdk1v(ji+1,jj,jk) & 341 !!$ +zdk1v(ji,jj,jk) + zdkv(ji+1,jj,jk) ) ) * fmask(ji,jj,jk) 342 !!$ ! j-flux at t-point 343 !!$ zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 344 !!$ 345 !!$ zmskt = 1._wp/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 346 !!$ + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp ) 347 !!$ 348 !!$ zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5_wp * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 349 !!$ 350 !!$ zjvt(ji,jj,jk) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & 351 !!$ + zcof2 * ( zdkv (ji,jj-1,jk) + zdk1v(ji,jj,jk) & 352 !!$ +zdk1v(ji,jj-1,jk) + zdkv (ji,jj,jk) ) ) * tmask(ji,jj,jk) 353 !!$ END DO 354 !!$ END DO 355 !!$ END DO 356 !!$ 357 !!$ ENDIF 358 359 360 361 IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) 362 363 DO jj = 2, jpjm1 364 DO ji = 2, jpim1 365 !DIR$ SHORTLOOP 366 ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the 367 ! If condition on jk>1. This does mean that there will be out-of-bounds *reads* 368 ! when jk is 1 but that doesn't matter. 369 !DIR$ SAFE_ADDRESS 370 DO jk = 1, jpkm1, 1 371 372 ! Vertical u- and v-shears at level jk and jk+1 373 ! --------------------------------------------- 374 ! surface boundary condition: zdku(jk=1)=zdku(jk=2) 375 ! zdkv(jk=1)=zdkv(jk=2) 376 zdk1u = ( ub(ji ,jj ,jk) -ub(ji ,jj ,jk+1) ) * umask(ji ,jj,jk+1) 377 zdk1i1u = ( ub(ji+1,jj ,jk) -ub(ji+1,jj ,jk+1) ) * umask(ji+1,jj,jk+1) 378 zdk1j1u = ( ub(ji ,jj+1,jk) -ub(ji ,jj+1,jk+1) ) * umask(ji ,jj+1,jk+1) 379 zdk1im1u = ( ub(ji-1,jj ,jk) -ub(ji-1,jj ,jk+1) ) * umask(ji-1,jj,jk+1) 380 zdk1jm1u = ( ub(ji ,jj-1,jk) -ub(ji ,jj-1,jk+1) ) * umask(ji ,jj-1,jk+1) 381 IF(jk > 1)THEN 382 zdku = ( ub(ji ,jj ,jk-1) - ub(ji ,jj ,jk) ) * umask(ji,jj,jk) 383 zdki1u = ( ub(ji+1,jj ,jk-1) - ub(ji+1,jj ,jk) ) * umask(ji+1,jj,jk) 384 zdkim1u= ( ub(ji-1,jj ,jk-1) - ub(ji-1,jj ,jk) ) * umask(ji-1,jj,jk) 385 zdkj1u = ( ub(ji ,jj+1,jk-1) - ub(ji ,jj+1,jk) ) * umask(ji,jj+1,jk) 386 zdkjm1u= ( ub(ji ,jj-1,jk-1) - ub(ji ,jj-1,jk) ) * umask(ji,jj-1,jk) 387 ELSE 388 zdku = zdk1u 389 zdki1u = zdk1i1u 390 zdkim1u= zdk1im1u 391 zdkj1u = zdk1j1u 392 zdkjm1u= zdk1jm1u 393 END IF 394 395 zdku = zdku + zdk1u 396 zdki1u = zdki1u + zdk1i1u 397 zdkim1u = zdkim1u + zdk1im1u 398 zdkj1u = zdkj1u +zdk1j1u 399 zdkjm1u = zdk1jm1u + zdkjm1u 400 401 ! volume elements 402 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 403 404 ! horizontal component of isopycnal momentum diffusive trends 405 406 ! -----f----- 407 ! Horizontal fluxes on U | 408 ! --------------------=== t u t 409 ! | 410 ! i-flux at t-point -----f----- 411 412 ! z-coordinate - partial steps : min(e3u) 413 aziuti1 = ( ((fsahmt(ji+1,jj,jk)+ahmb0) * e2t(ji+1,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji+1,jj)) * & 414 ( ub(ji+1,jj,jk) - ub(ji,jj,jk) ) + (- aht0 * e2t(ji+1,jj) * & 415 (1._wp/MAX( umask(ji,jj,jk )+umask(ji+1,jj,jk+1) & 416 + umask(ji,jj,jk+1)+umask(ji+1,jj,jk ), 1._wp )) * 0.5 * & 417 ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) )) * & 418 ( zdku + zdki1u ) ) * tmask(ji+1,jj,jk) 419 420 aziut = ( ((fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * & 421 ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) + (- aht0 * e2t(ji,jj) * & 422 (1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 423 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * & 424 ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )) * & 425 ( zdku + zdkim1u ) ) * tmask(ji,jj,jk) 426 427 ! j-flux at f-point - identical to non-ln_zps 428 429 azjuf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)) * & 430 ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) + (- aht0 * e1f(ji,jj) * & 431 (1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & 432 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * & 433 ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )) * & 434 ( zdku + zdkj1u ) ) * fmask(ji,jj,jk) 435 436 azjufjm1 = ( (( fsahmf(ji,jj-1,jk) + ahmb0 ) * e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1)) * & 437 ( ub(ji,jj,jk) - ub(ji,jj-1,jk) ) + (- aht0 * e1f(ji,jj-1) * & 438 (1./MAX( umask(ji,jj,jk )+umask(ji,jj-1,jk+1) & 439 + umask(ji,jj,jk+1)+umask(ji,jj-1,jk ), 1. )) * 0.5 * & 440 ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) )) * & 441 ( zdku + zdkjm1u ) ) * fmask(ji,jj-1,jk) 442 443 ! Second derivative (divergence) and add to the general trend 444 ! ----------------------------------------------------------- 445 446 zuah =( & 447 ! ziut (ji+1,jj,jk) - & 448 aziuti1 - & 449 ! ziut (ji ,jj,jk) + & 450 aziut + & 451 ! zjuf (ji ,jj,jk) - & 452 azjuf - & 453 ! zjuf (ji,jj-1,jk) & 454 azjufjm1 & 455 ) / zbu 456 457 ! add the trends to the general trends 458 ua (ji,jj,jk) = ua (ji,jj,jk) + zuah 459 460 END DO 461 462 !DIR$ SHORTLOOP 463 ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the 464 ! If condition on jk>1. This does mean that there will be out-of-bounds *reads* 465 ! when jk is 1 but that doesn't matter. 466 !DIR$ SAFE_ADDRESS 467 DO jk = 1, jpkm1, 1 468 469 zdk1v = ( vb(ji ,jj ,jk) -vb(ji ,jj ,jk+1) ) * vmask(ji ,jj,jk+1) 470 zdk1i1v = ( vb(ji+1,jj ,jk) -vb(ji+1,jj ,jk+1) ) * vmask(ji+1,jj,jk+1) 471 zdk1im1v = ( vb(ji-1,jj ,jk) -vb(ji-1,jj ,jk+1) ) * vmask(ji-1,jj,jk+1) 472 zdk1j1v = ( vb(ji ,jj+1,jk) -vb(ji ,jj+1,jk+1) ) * vmask(ji ,jj+1,jk+1) 473 zdk1jm1v = ( vb(ji ,jj-1,jk) -vb(ji ,jj-1,jk+1) ) * vmask(ji ,jj-1,jk+1) 474 IF(jk > 1)THEN 475 zdkv = ( vb(ji ,jj ,jk-1) - vb(ji ,jj ,jk) ) * vmask(ji,jj,jk) 476 zdki1v = ( vb(ji+1,jj ,jk-1) - vb(ji+1,jj ,jk) ) * vmask(ji+1,jj,jk) 477 zdkim1v= ( vb(ji-1,jj ,jk-1) - vb(ji-1,jj ,jk) ) * vmask(ji-1,jj,jk) 478 zdkj1v = ( vb(ji ,jj+1,jk-1) - vb(ji ,jj+1,jk) ) * vmask(ji,jj+1,jk) 479 zdkjm1v= ( vb(ji ,jj-1,jk-1) - vb(ji ,jj-1,jk) ) * vmask(ji,jj-1,jk) 480 ELSE 481 zdkv = zdk1v 482 zdki1v = zdk1i1v 483 zdkim1v= zdk1im1v 484 zdkj1v = zdk1j1v 485 zdkjm1v= zdk1jm1v 486 END IF 487 488 zdkv = zdkv + zdk1v 489 zdki1v = zdk1i1v + zdki1v 490 zdkim1v = zdkim1v +zdk1im1v 491 zdkj1v = zdk1j1v + zdkj1v 492 zdkjm1v = zdkjm1v +zdk1jm1v 493 494 ! | t | 495 ! Horizontal fluxes on V | | 496 ! --------------------=== f---v---f 497 ! | | 498 ! | t | 499 500 ! i-flux at f-point - identical to non-ln_zps case 501 azivf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / & 502 e1f(ji,jj)) * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & 503 + (- aht0 * e2f(ji,jj) / & 504 MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 505 + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp ) * & 506 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * & 507 ( zdkv + zdki1v ) ) * fmask(ji,jj,jk) 508 509 azivfim1 = ( (( fsahmf(ji-1,jj,jk) + ahmb0 ) * e2f(ji-1,jj) * & 510 fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( vb(ji,jj,jk) - vb(ji-1,jj,jk) ) & 511 + (- aht0 * e2f(ji-1,jj) / & 512 MAX( vmask(ji,jj,jk )+vmask(ji-1,jj,jk+1) & 513 + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk ), 1._wp ) * 0.5_wp * & 514 ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * & 515 ( zdkim1v + zdkv ) ) * fmask(ji-1,jj,jk) 516 517 ! j-flux at t-point - min(e3u) instead of e3t 518 azjvtj1 = ( & 519 ((fsahmt(ji,jj+1,jk)+ahmb0) * e1t(ji,jj+1) * MIN( fse3v(ji,jj+1,jk), fse3v(ji,jj,jk) ) / & 520 e2t(ji,jj+1)) * ( vb(ji,jj+1,jk) - vb(ji,jj,jk) ) & 521 + (- aht0 * e1t(ji,jj+1) / & 522 MAX( vmask(ji,jj,jk )+vmask(ji,jj+1,jk+1) & 523 + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk ), 1._wp ) * & 524 0.5_wp * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * & 525 ( zdkv + zdkj1v ) & 526 ) * tmask(ji,jj+1,jk) 527 528 azjvt = ( & 529 ((fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / & 530 e2t(ji,jj)) * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & 531 + (- aht0 * e1t(ji,jj) /MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 532 + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp ) * 0.5_wp * & 533 ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * & 534 ( zdkjm1v + zdkv ) ) * tmask(ji,jj,jk) 535 536 ! volume elements 537 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 538 539 540 ! Second derivative (divergence) and add to the general trend 541 ! ----------------------------------------------------------- 542 zvah =( & 543 ! zivf (ji ,jj ,jk) - & 544 azivf - & 545 ! zivf (ji-1,jj ,jk) + & 546 azivfim1 + & 547 ! zjvt (ji ,jj+1,jk) - & 548 azjvtj1 - & 549 ! zjvt (ji ,jj ,jk) & 550 azjvt & 551 ) / zbv 552 553 ! add the trend to the general trend 554 va (ji,jj,jk) = va (ji,jj,jk) + zvah 555 END DO 556 END DO 557 END DO 558 559 ELSE ! other coordinate system (zco or sco) : e3t 560 561 DO jj = 2, jpjm1 562 DO ji = 2, jpim1 563 !DIR$ SHORTLOOP 564 ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the 565 ! If condition on jk>1. This does mean that there will be out-of-bounds *reads* 566 ! when jk is 1 but that doesn't matter. 567 !DIR$ SAFE_ADDRESS 568 DO jk = 1, jpkm1, 1 569 570 ! Vertical u- and v-shears at level jk and jk+1 571 ! --------------------------------------------- 572 ! surface boundary condition: zdku(jk=1)=zdku(jk=2) 573 ! zdkv(jk=1)=zdkv(jk=2) 574 zdk1u = ( ub(ji ,jj ,jk) -ub(ji ,jj ,jk+1) ) * umask(ji ,jj,jk+1) 575 zdk1i1u = ( ub(ji+1,jj ,jk) -ub(ji+1,jj ,jk+1) ) * umask(ji+1,jj,jk+1) 576 zdk1j1u = ( ub(ji ,jj+1,jk) -ub(ji ,jj+1,jk+1) ) * umask(ji ,jj+1,jk+1) 577 zdk1im1u = ( ub(ji-1,jj ,jk) -ub(ji-1,jj ,jk+1) ) * umask(ji-1,jj,jk+1) 578 zdk1jm1u = ( ub(ji ,jj-1,jk) -ub(ji ,jj-1,jk+1) ) * umask(ji ,jj-1,jk+1) 579 IF(jk > 1)THEN 580 zdku = ( ub(ji ,jj ,jk-1) - ub(ji ,jj ,jk) ) * umask(ji,jj,jk) 581 zdki1u = ( ub(ji+1,jj ,jk-1) - ub(ji+1,jj ,jk) ) * umask(ji+1,jj,jk) 582 zdkim1u= ( ub(ji-1,jj ,jk-1) - ub(ji-1,jj ,jk) ) * umask(ji-1,jj,jk) 583 zdkj1u = ( ub(ji ,jj+1,jk-1) - ub(ji ,jj+1,jk) ) * umask(ji,jj+1,jk) 584 zdkjm1u= ( ub(ji ,jj-1,jk-1) - ub(ji ,jj-1,jk) ) * umask(ji,jj-1,jk) 585 ELSE 586 zdku = zdk1u 587 zdki1u = zdk1i1u 588 zdkim1u= zdk1im1u 589 zdkj1u = zdk1j1u 590 zdkjm1u= zdk1jm1u 591 END IF 592 593 zdku = zdku + zdk1u 594 zdki1u = zdki1u + zdk1i1u 595 zdkim1u = zdkim1u + zdk1im1u 596 zdkj1u = zdkj1u +zdk1j1u 597 zdkjm1u = zdk1jm1u + zdkjm1u 598 599 ! volume element 600 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 601 602 ! horizontal component of isopycnal momentum diffusive trends 603 604 ! -----f----- 605 ! Horizontal fluxes on U | 606 ! --------------------=== t u t 607 ! | 608 ! i-flux at t-point -----f----- 609 610 ! other coordinate system (zco or sco) : e3t 611 aziuti1 = ( ((fsahmt(ji+1,jj,jk)+ahmb0) * e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * & 612 ( ub(ji+1,jj,jk) - ub(ji,jj,jk) ) + (- aht0 * e2t(ji+1,jj) * & 613 (1./MAX( umask(ji,jj,jk )+umask(ji+1,jj,jk+1) & 614 + umask(ji,jj,jk+1)+umask(ji+1,jj,jk ), 1. )) * 0.5 * & 615 ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) )) * & 616 ( zdku + zdki1u ) ) * tmask(ji+1,jj,jk) 617 618 aziut = ( ((fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * & 619 ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) + (- aht0 * e2t(ji,jj) * & 620 (1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 621 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * & 622 ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )) * & 623 ( zdku + zdkim1u ) ) * tmask(ji,jj,jk) 624 625 azjuf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)) * & 626 ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) + (- aht0 * e1f(ji,jj) * & 627 (1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & 628 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * & 629 ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )) * & 630 ( zdku + zdkj1u ) ) * fmask(ji,jj,jk) 631 632 azjufjm1 = ( (( fsahmf(ji,jj-1,jk) + ahmb0 ) * e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1)) * & 633 ( ub(ji,jj,jk) - ub(ji,jj-1,jk) ) + (- aht0 * e1f(ji,jj-1) * & 634 (1./MAX( umask(ji,jj,jk )+umask(ji,jj-1,jk+1) & 635 + umask(ji,jj,jk+1)+umask(ji,jj-1,jk ), 1. )) * 0.5 * & 636 ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) )) * & 637 ( zdku + zdkjm1u ) ) * fmask(ji,jj-1,jk) 638 639 640 ! Second derivative (divergence) and add to the general trend 641 ! ----------------------------------------------------------- 642 zuah =( & 643 ! ziut (ji+1,jj,jk) - & 644 aziuti1 - & 645 ! ziut (ji ,jj,jk) + & 646 aziut + & 647 ! zjuf (ji ,jj,jk) - & 648 azjuf - & 649 ! zjuf (ji,jj-1,jk) & 650 azjufjm1 & 651 ) / zbu 652 653 ! add the trend to the general trend 654 ua (ji,jj,jk) = ua (ji,jj,jk) + zuah 655 656 END DO 657 658 !DIR$ SHORTLOOP 659 ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the 660 ! If condition on jk>1. This does mean that there will be out-of-bounds *reads* 661 ! when jk is 1 but that doesn't matter. 662 !DIR$ SAFE_ADDRESS 663 DO jk = 1, jpkm1, 1 664 665 zdk1v = ( vb(ji ,jj ,jk) -vb(ji ,jj ,jk+1) ) * vmask(ji ,jj,jk+1) 666 zdk1i1v = ( vb(ji+1,jj ,jk) -vb(ji+1,jj ,jk+1) ) * vmask(ji+1,jj,jk+1) 667 zdk1im1v = ( vb(ji-1,jj ,jk) -vb(ji-1,jj ,jk+1) ) * vmask(ji-1,jj,jk+1) 668 zdk1j1v = ( vb(ji ,jj+1,jk) -vb(ji ,jj+1,jk+1) ) * vmask(ji ,jj+1,jk+1) 669 zdk1jm1v = ( vb(ji ,jj-1,jk) -vb(ji ,jj-1,jk+1) ) * vmask(ji ,jj-1,jk+1) 670 IF(jk > 1)THEN 671 zdkv = ( vb(ji ,jj ,jk-1) - vb(ji ,jj ,jk) ) * vmask(ji,jj,jk) 672 zdki1v = ( vb(ji+1,jj ,jk-1) - vb(ji+1,jj ,jk) ) * vmask(ji+1,jj,jk) 673 zdkim1v= ( vb(ji-1,jj ,jk-1) - vb(ji-1,jj ,jk) ) * vmask(ji-1,jj,jk) 674 zdkj1v = ( vb(ji ,jj+1,jk-1) - vb(ji ,jj+1,jk) ) * vmask(ji,jj+1,jk) 675 zdkjm1v= ( vb(ji ,jj-1,jk-1) - vb(ji ,jj-1,jk) ) * vmask(ji,jj-1,jk) 676 ELSE 677 zdkv = zdk1v 678 zdki1v = zdk1i1v 679 zdkim1v= zdk1im1v 680 zdkj1v = zdk1j1v 681 zdkjm1v= zdk1jm1v 682 END IF 683 684 zdkv = zdkv + zdk1v 685 zdki1v = zdk1i1v + zdki1v 686 zdkim1v = zdkim1v +zdk1im1v 687 zdkj1v = zdk1j1v + zdkj1v 688 zdkjm1v = zdkjm1v +zdk1jm1v 689 690 azivf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / & 691 e1f(ji,jj)) * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & 692 + (- aht0 * e2f(ji,jj) / & 693 MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 694 + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp ) * & 695 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * & 696 ( zdkv + zdki1v ) ) * fmask(ji,jj,jk) 697 698 azivfim1 = ( (( fsahmf(ji-1,jj,jk) + ahmb0 ) * e2f(ji-1,jj) * & 699 fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( vb(ji,jj,jk) - vb(ji-1,jj,jk) ) & 700 + (- aht0 * e2f(ji-1,jj) / & 701 MAX( vmask(ji,jj,jk )+vmask(ji-1,jj,jk+1) & 702 + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk ), 1._wp ) * 0.5_wp * & 703 ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * & 704 ( zdkim1v + zdkv ) ) * fmask(ji-1,jj,jk) 705 706 azjvtj1 = ( & 707 ((fsahmt(ji,jj+1,jk)+ahmb0) * e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / & 708 e2t(ji,jj+1)) * ( vb(ji,jj+1,jk) - vb(ji,jj,jk) ) & 709 + (- aht0 * e1t(ji,jj+1) / & 710 MAX( vmask(ji,jj,jk )+vmask(ji,jj+1,jk+1) & 711 + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk ), 1._wp ) * & 712 0.5_wp * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * & 713 ( zdkv + zdkj1v ) & 714 ) * tmask(ji,jj+1,jk) 715 716 azjvt = ( ((fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * & 717 ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & 718 + (- aht0 * e1t(ji,jj) /MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 719 + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp ) * 0.5_wp * & 720 ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * & 721 ( zdkjm1v + zdkv ) ) * tmask(ji,jj,jk) 722 723 ! volume elements 724 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 725 726 ! Second derivative (divergence) and add to the general trend 727 ! ----------------------------------------------------------- 728 zvah =( & 729 ! zivf (ji ,jj ,jk) - & 730 azivf - & 731 ! zivf (ji-1,jj ,jk) + & 732 azivfim1 + & 733 ! zjvt (ji ,jj+1,jk) - & 734 azjvtj1 - & 735 ! zjvt (ji ,jj ,jk) & 736 azjvt & 737 ) / zbv 738 739 ! add the trend to the general trend 740 va (ji,jj,jk) = va (ji,jj,jk) + zvah 741 END DO 742 END DO 743 END DO 744 745 END IF ! ln_zps 746 #else 172 747 ! ! =============== 173 748 DO jk = 1, jpkm1 ! Horizontal slab … … 320 895 END DO ! End of slab 321 896 ! ! =============== 897 #endif 898 CALL timing_stop('dyn_ldf_iso_hslab','section') 322 899 323 900 ! print sum trends (used for debugging) … … 325 902 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 326 903 327 904 CALL timing_start('dyn_ldf_iso_vslab') 905 906 #if defined key_z_first 907 ! ! =============== 908 DO jj = 2, jpjm1 ! Vertical slab 909 ! ! =============== 910 911 ! I. vertical trends associated with the lateral mixing 912 ! ===================================================== 913 ! (excluding the vertical flux proportional to dk[t] 914 915 916 ! I.1 horizontal momentum gradient 917 ! -------------------------------- 918 919 DO ji = 2, jpi 920 DO jk = 1, jpk 921 ! i-gradient of u at jj 922 zdiu (ji,jk) = tmask(ji,jj ,jk) * ( ub(ji,jj ,jk) - ub(ji-1,jj ,jk) ) 923 ! j-gradient of u and v at jj 924 zdju (ji,jk) = fmask(ji,jj ,jk) * ( ub(ji,jj+1,jk) - ub(ji ,jj ,jk) ) 925 zdjv (ji,jk) = tmask(ji,jj ,jk) * ( vb(ji,jj ,jk) - vb(ji ,jj-1,jk) ) 926 ! j-gradient of u and v at jj+1 927 zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj ,jk) - ub(ji ,jj-1,jk) ) 928 zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji ,jj ,jk) ) 929 END DO 930 END DO 931 932 DO ji = 1, jpim1 933 DO jk = 1, jpk 934 ! i-gradient of v at jj 935 zdiv (ji,jk) = fmask(ji,jj ,jk) * ( vb(ji+1,jj,jk) - vb(ji ,jj ,jk) ) 936 END DO 937 END DO 938 939 940 ! I.2 Vertical fluxes 941 ! ------------------- 942 943 ! Surface and bottom vertical fluxes set to zero 944 !!$ DO ji = 1, jpi 945 !!$ zfuw(ji, 1 ) = 0.e0 946 !!$ zfvw(ji, 1 ) = 0.e0 947 !!$ zfuw(ji,jpk) = 0.e0 948 !!$ zfvw(ji,jpk) = 0.e0 949 !!$ END DO 950 951 ! interior (2=<jk=<jpk-1) on U field 952 DO ji = 2, jpim1 953 DO jk = 2, jpkm1 954 zcoef0= 0.5 * aht0 * umask(ji,jj,jk) 955 956 zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 957 zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) ) 958 959 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) & 960 + tmask(ji,jj,jk )+tmask(ji+1,jj,jk ), 1. ) 961 zmkf = 1./MAX( fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1) & 962 + fmask(ji,jj-1,jk )+fmask(ji,jj,jk ), 1. ) 963 964 zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi 965 zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj 966 ! vertical flux on u field 967 zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) & 968 +zdiu (ji,jk ) + zdiu (ji+1,jk ) ) & 969 + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) & 970 +zdj1u(ji,jk ) + zdju (ji ,jk ) ) 971 ! update avmu (add isopycnal vertical coefficient to avmu) 972 ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 973 avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0 974 END DO 975 END DO 976 977 ! interior (2=<jk=<jpk-1) on V field 978 DO ji = 2, jpim1 979 DO jk = 2, jpkm1 980 zcoef0= 0.5 * aht0 * vmask(ji,jj,jk) 981 982 zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) 983 zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) ) 984 985 zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) & 986 + fmask(ji-1,jj,jk )+fmask(ji,jj,jk ), 1. ) 987 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) & 988 + tmask(ji,jj,jk )+tmask(ji,jj+1,jk ), 1. ) 989 990 zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi 991 zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj 992 ! vertical flux on v field 993 zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) & 994 +zdiv (ji,jk ) + zdiv (ji-1,jk ) ) & 995 + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) & 996 +zdjv (ji,jk ) + zdj1v(ji ,jk ) ) 997 ! update avmv (add isopycnal vertical coefficient to avmv) 998 ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 999 avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0 1000 END DO 1001 ! END DO 1002 1003 1004 ! I.3 Divergence of vertical fluxes added to the general tracer trend 1005 ! ------------------------------------------------------------------- 1006 ! DO ji = 2, jpim1 1007 zfuw(ji,1) = 0.0_wp 1008 zfuw(ji,jpk) = 0.0_wp 1009 zfvw(ji,1) = 0.0_wp 1010 zfvw(ji,jpk) = 0.0_wp 1011 1012 DO jk = 1, jpkm1 1013 ! volume elements 1014 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 1015 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 1016 ! part of the k-component of isopycnal momentum diffusive trends 1017 zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu 1018 zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv 1019 ! add the trends to the general trends 1020 ua(ji,jj,jk) = ua(ji,jj,jk) + zuav 1021 va(ji,jj,jk) = va(ji,jj,jk) + zvav 1022 END DO 1023 END DO 1024 ! ! =============== 1025 END DO ! End of vertical slab 1026 ! ! =============== 1027 #else 328 1028 ! ! =============== 329 1029 DO jj = 2, jpjm1 ! Vertical slab … … 441 1141 END DO ! End of slab 442 1142 ! ! =============== 1143 #endif 1144 CALL timing_stop('dyn_ldf_iso_vslab','section') 443 1145 444 1146 IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) ) CALL ctl_stop('dyn_ldf_iso: failed to release workspace arrays') 1147 ! 1148 CALL timing_stop('dyn_ldf_iso','section') 445 1149 ! 446 1150 END SUBROUTINE dyn_ldf_iso -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r3211 r3432 281 281 ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 282 282 zs_t (:,:) = e1t(:,:) * e2t(:,:) 283 zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)284 zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)283 zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) ) 284 zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) ) 285 285 ! 286 286 IF( ln_dynadv_vec ) THEN -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r3211 r3432 122 122 USE wrk_nemo, ONLY: zun => wrk_2d_14 , zun_e => wrk_2d_15 , zub_e => wrk_2d_16 , zu_sum => wrk_2d_17 123 123 USE wrk_nemo, ONLY: zvn => wrk_2d_18 , zvn_e => wrk_2d_19 , zvb_e => wrk_2d_20 , zv_sum => wrk_2d_21 124 USE timing, ONLY: timing_start, timing_stop 124 125 ! 125 126 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 132 133 REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - 133 134 REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - 135 ! DCSE_NEMO - this temporary array used to take a copy of the surface values 136 ! of the rhd array. 137 REAL(wp), ALLOCATABLE :: zrhd_1(:,:) ! local array to hold surface 138 ! values of rhd in z-first case 134 139 !!---------------------------------------------------------------------- 140 141 CALL timing_start('dyn_spg_ts') 135 142 136 143 IF( wrk_in_use(2, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, & … … 178 185 zv_sld = 0.e0 ; zv_asp = 0.e0 179 186 187 #if defined key_z_first 188 ! DCSE_NEMO - allocate temp. array and then use it to store surface values 189 ! (jk=1) of full rhd array. 190 IF( lk_vvl) THEN 191 ALLOCATE(zrhd_1(jpi,jpj)) 192 DO jj=1,jpj,1 193 DO ji=1,jpi,1 194 zrhd_1(ji,jj) = rhd(ji,jj,1) 195 END DO 196 END DO 197 END IF 198 #endif 180 199 ! ----------------------------------------------------------------------------- 181 200 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) … … 280 299 DO jj = 2, jpjm1 281 300 DO ji = fs_2, fs_jpim1 ! vector opt. 301 #if defined key_z_first 302 zcu(ji,jj) = zcu(ji,jj) - grav * ( ( zrhd_1(ji+1,jj ) + 1 ) * sshn(ji+1,jj ) & 303 & - ( zrhd_1(ji ,jj ) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 304 zcv(ji,jj) = zcv(ji,jj) - grav * ( ( zrhd_1(ji ,jj+1) + 1 ) * sshn(ji ,jj+1) & 305 & - ( zrhd_1(ji ,jj ) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 306 #else 282 307 zcu(ji,jj) = zcu(ji,jj) - grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & 283 308 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 284 309 zcv(ji,jj) = zcv(ji,jj) - grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 285 310 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 311 #endif 286 312 END DO 287 313 END DO … … 317 343 DO jj = 2, jpjm1 318 344 DO ji = fs_2, fs_jpim1 ! vector opt. 345 #if defined key_z_first 346 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) & 347 & / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.0_wp - umask_1(ji,jj) ) 348 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) & 349 & / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.0_wp - vmask_1(ji,jj) ) 350 #else 319 351 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) & 320 352 & / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) ) 321 353 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) & 322 354 & / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) ) 355 #endif 323 356 END DO 324 357 END DO … … 338 371 ! 339 372 IF( lk_vvl ) THEN 373 #if defined key_z_first 374 ub_b(:,:) = ub_b(:,:) * umask_1(:,:) / ( hu_0(:,:) + sshu_b(:,:) + 1.0_wp - umask_1(:,:) ) 375 vb_b(:,:) = vb_b(:,:) * vmask_1(:,:) / ( hv_0(:,:) + sshv_b(:,:) + 1.0_wp - vmask_1(:,:) ) 376 #else 340 377 ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) ) 341 378 vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) ) 379 #endif 342 380 ELSE 343 381 ub_b(:,:) = ub_b(:,:) * hur(:,:) … … 417 455 DO jj = 2, jpjm1 ! leap-frog on ssh_e 418 456 DO ji = fs_2, fs_jpim1 ! vector opt. 457 #if defined key_z_first 458 ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask_1(ji,jj) 459 #else 419 460 ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 461 #endif 420 462 END DO 421 463 END DO … … 431 473 ! surface pressure gradient 432 474 IF( lk_vvl) THEN 475 #if defined key_z_first 476 zu_spg = -grav * ( ( zrhd_1(ji+1,jj ) + 1 ) * sshn_e(ji+1,jj ) & 477 & - ( zrhd_1(ji ,jj ) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 478 zv_spg = -grav * ( ( zrhd_1(ji ,jj+1) + 1 ) * sshn_e(ji ,jj+1) & 479 & - ( zrhd_1(ji ,jj ) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 480 #else 433 481 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 434 482 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 435 483 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 436 484 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 485 #endif 437 486 ELSE 438 487 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) … … 447 496 zv_cor =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 448 497 ! after velocities with implicit bottom friction 498 #if defined key_z_first 499 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask_1(ji,jj) & 500 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 501 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask_1(ji,jj) & 502 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 503 #else 449 504 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 450 505 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 451 506 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 452 507 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 508 #endif 453 509 END DO 454 510 END DO … … 459 515 ! surface pressure gradient 460 516 IF( lk_vvl) THEN 517 #if defined key_z_first 518 zu_spg = -grav * ( ( zrhd_1(ji+1,jj ) + 1 ) * sshn_e(ji+1,jj ) & 519 & - ( zrhd_1(ji ,jj ) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 520 zv_spg = -grav * ( ( zrhd_1(ji ,jj+1) + 1 ) * sshn_e(ji ,jj+1) & 521 & - ( zrhd_1(ji ,jj ) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 522 #else 461 523 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 462 524 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 463 525 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 464 526 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 527 #endif 465 528 ELSE 466 529 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) … … 473 536 zv_cor = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj) 474 537 ! after velocities with implicit bottom friction 538 #if defined key_z_first 539 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask_1(ji,jj) & 540 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 541 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask_1(ji,jj) & 542 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 543 #else 475 544 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 476 545 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 477 546 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 478 547 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 548 #endif 479 549 END DO 480 550 END DO … … 485 555 ! surface pressure gradient 486 556 IF( lk_vvl) THEN 557 #if defined key_z_first 558 zu_spg = -grav * ( ( zrhd_1(ji+1,jj ) + 1 ) * sshn_e(ji+1,jj ) & 559 & - ( zrhd_1(ji ,jj ) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 560 zv_spg = -grav * ( ( zrhd_1(ji ,jj+1) + 1 ) * sshn_e(ji ,jj+1) & 561 & - ( zrhd_1(ji ,jj ) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 562 #else 487 563 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 488 564 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 489 565 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 490 566 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 567 #endif 491 568 ELSE 492 569 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) … … 499 576 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj) 500 577 ! after velocities with implicit bottom friction 578 #if defined key_z_first 579 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask_1(ji,jj) & 580 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 581 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask_1(ji,jj) & 582 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 583 #else 501 584 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 502 585 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 503 586 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 504 587 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 588 #endif 505 589 END DO 506 590 END DO … … 546 630 DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points 547 631 DO ji = 1, fs_jpim1 ! Vector opt. 632 #if defined key_z_first 633 zsshun_e(ji,jj) = 0.5 * umask_1(ji,jj) / ( e1u(ji,jj) * e2u(ji,jj) ) & 634 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & 635 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 636 zsshvn_e(ji,jj) = 0.5 * vmask_1(ji,jj) / ( e1v(ji,jj) * e2v(ji,jj) ) & 637 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 638 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 639 #else 548 640 zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 549 641 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & … … 552 644 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 553 645 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 646 #endif 554 647 END DO 555 648 END DO … … 559 652 hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points 560 653 hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 654 #if defined key_z_first 655 DO jj = 1,jpj,1 656 DO ji=1,jpi,1 657 hur_e(ji,jj) = umask_1(ji,jj) / ( hu_e(ji,jj) + 1.0_wp - umask_1(ji,jj) ) 658 hvr_e(ji,jj) = vmask_1(ji,jj) / ( hv_e(ji,jj) + 1.0_wp - vmask_1(ji,jj) ) 659 END DO 660 END DO 661 #else 561 662 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) ) 562 663 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) ) 664 #endif 563 665 ! 564 666 ENDIF … … 584 686 ! 585 687 ! !* update the general momentum trend 688 #if defined key_z_first 689 DO jj=1,jpj,1 690 DO ji=1,jpi,1 691 zu_asp = ( zu_sum(ji,jj) - ub_b(ji,jj) ) / z2dt_b 692 zv_asp = ( zv_sum(ji,jj) - vb_b(ji,jj) ) / z2dt_b 693 DO jk=1,jpkm1 694 ua(ji,jj,jk) = ua(ji,jj,jk) + zu_asp 695 va(ji,jj,jk) = va(ji,jj,jk) + zv_asp 696 END DO 697 END DO 698 END DO 699 #else 586 700 DO jk=1,jpkm1 587 701 ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b 588 702 va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b 589 703 END DO 704 #endif 590 705 un_b (:,:) = zu_sum(:,:) 591 706 vn_b (:,:) = zv_sum(:,:) … … 599 714 & 11,12,13,14,15,16,17,18,19,20,21) ) & 600 715 CALL ctl_stop('dyn_spg_ts: failed to release workspace arrays') 716 ! 717 CALL timing_stop('dyn_spg_ts','section') 601 718 ! 602 719 END SUBROUTINE dyn_spg_ts -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r3211 r3432 157 157 ! !------------------------------------------! 158 158 #if defined key_z_first 159 fsdept(:,:,1:jpkm1) = fsdept_n(:,:,1:jpkm1) ! now local depths stored in fsdep. arrays 160 fsdepw(:,:,1:jpkm1) = fsdepw_n(:,:,1:jpkm1) 161 fsde3w(:,:,1:jpkm1) = fsde3w_n(:,:,1:jpkm1) 162 ! 163 fse3t (:,:,1:jpkm1) = fse3t_n (:,:,1:jpkm1) ! vertical scale factors stored in fse3. arrays 164 fse3u (:,:,1:jpkm1) = fse3u_n (:,:,1:jpkm1) 165 fse3v (:,:,1:jpkm1) = fse3v_n (:,:,1:jpkm1) 166 fse3f (:,:,1:jpkm1) = fse3f_n (:,:,1:jpkm1) 167 fse3w (:,:,1:jpkm1) = fse3w_n (:,:,1:jpkm1) 168 fse3uw(:,:,1:jpkm1) = fse3uw_n(:,:,1:jpkm1) 169 fse3vw(:,:,1:jpkm1) = fse3vw_n(:,:,1:jpkm1) 159 ! DCSE_NEMO: can't use implicit loop over k here because the domzgr_substitute.h90 160 ! file causes the line below to be expanded to: 161 ! gdept_1(1:jpkm1,:,:) = (gdept(1:jpkm1,:,:)*(1.+sshn(:,:)*mut(1:jpkm1,:,:))) 162 ! which contains non-conforming array expressions. 163 DO jj=1,jpj,1 164 DO ji=1,jpi,1 165 DO jk=1,jpk,1 166 fsdept(ji,jj,jk) = fsdept_n(ji,jj,jk) ! now local depths stored in fsdep. arrays 167 END DO 168 END DO 169 END DO 170 DO jj=1,jpj,1 171 DO ji=1,jpi,1 172 DO jk=1,jpk,1 173 fsdepw(ji,jj,jk) = fsdepw_n(ji,jj,jk) 174 END DO 175 END DO 176 END DO 177 DO jj=1,jpj,1 178 DO ji=1,jpi,1 179 DO jk=1,jpk,1 180 fsde3w(ji,jj,jk) = fsde3w_n(ji,jj,jk) 181 END DO 182 END DO 183 END DO 184 ! 185 DO jj=1,jpj,1 186 DO ji=1,jpi,1 187 DO jk=1,jpk,1 188 fse3t (ji,jj,jk) = fse3t_n (ji,jj,jk) ! vertical scale factors stored in fse3. arrays 189 END DO 190 END DO 191 END DO 192 DO jj=1,jpj,1 193 DO ji=1,jpi,1 194 DO jk=1,jpk,1 195 fse3u (ji,jj,jk) = fse3u_n (ji,jj,jk) 196 END DO 197 END DO 198 END DO 199 DO jj=1,jpj,1 200 DO ji=1,jpi,1 201 DO jk=1,jpk,1 202 fse3v (ji,jj,jk) = fse3v_n (ji,jj,jk) 203 END DO 204 END DO 205 END DO 206 DO jj=1,jpj,1 207 DO ji=1,jpi,1 208 DO jk=1,jpk,1 209 fse3f (ji,jj,jk) = fse3f_n (ji,jj,jk) 210 END DO 211 END DO 212 END DO 213 DO jj=1,jpj,1 214 DO ji=1,jpi,1 215 DO jk=1,jpk,1 216 fse3w (ji,jj,jk) = fse3w_n (ji,jj,jk) 217 END DO 218 END DO 219 END DO 220 221 222 DO jj=1,jpj,1 223 DO ji=1,jpi,1 224 DO jk=1,jpk,1 225 fse3uw(ji,jj,jk) = fse3uw_n(ji,jj,jk) 226 END DO 227 END DO 228 END DO 229 230 DO jj=1,jpj,1 231 DO ji=1,jpi,1 232 DO jk=1,jpk,1 233 fse3vw(ji,jj,jk) = fse3vw_n(ji,jj,jk) 234 END DO 235 END DO 236 END DO 170 237 #else 171 238 DO jk = 1, jpkm1 -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
r3211 r3432 35 35 # endif 36 36 USE zpermute, ONLY : permute_z_last ! Re-order a 3d array back to external (z-last) ordering 37 37 USE timing, ONLY : timing_start, timing_stop 38 38 IMPLICIT NONE 39 39 PUBLIC ! must be public to be able to access iom_def through iom … … 184 184 ! end halo size for x,y dimensions 185 185 !--------------------------------------------------------------------- 186 187 CALL timing_start('iom_open') 188 186 189 ! Initializations and control 187 190 ! ============= … … 323 326 ENDIF 324 327 ! 328 CALL timing_stop('iom_open') 329 325 330 END SUBROUTINE iom_open 326 331 … … 340 345 CHARACTER(LEN=100) :: clinfo ! info character 341 346 !--------------------------------------------------------------------- 347 ! 348 CALL timing_start('iom_close') 342 349 ! 343 350 clinfo = ' iom_close ~~~ ' … … 373 380 ENDIF 374 381 ! 382 CALL timing_stop('iom_close') 383 375 384 END SUBROUTINE iom_close 376 385 … … 392 401 LOGICAL :: llstop ! local definition of ldstop 393 402 !!----------------------------------------------------------------------- 403 CALL timing_start('iom_varid') 404 394 405 iom_varid = 0 ! default definition 395 406 ! do we call ctl_stop if we look for non-existing variable? … … 441 452 ENDIF 442 453 ! 454 CALL timing_stop('iom_varid') 455 ! 443 456 END FUNCTION iom_varid 444 457 … … 455 468 ! 456 469 IF( kiomid > 0 ) THEN 470 CALL timing_start('iom_g0d') 457 471 idvar = iom_varid( kiomid, cdvar ) 458 472 IF( iom_file(kiomid)%nfid > 0 .AND. idvar > 0 ) THEN … … 465 479 END SELECT 466 480 ENDIF 481 CALL timing_stop('iom_g0d') 467 482 ENDIF 468 483 END SUBROUTINE iom_g0d 469 484 470 SUBROUTINE iom_g1d( kiomid, kdom, cdvar, pvar, ktime, kstart, kcount )485 SUBROUTINE iom_g1d( kiomid, kdom, cdvar, pvar, ktime, kstart, kcount, lzfirst ) 471 486 INTEGER , INTENT(in ) :: kiomid ! Identifier of the file 472 487 INTEGER , INTENT(in ) :: kdom ! Type of domain to be read … … 476 491 INTEGER , INTENT(in ), DIMENSION(1), OPTIONAL :: kstart ! start axis position of the reading 477 492 INTEGER , INTENT(in ), DIMENSION(1), OPTIONAL :: kcount ! number of points in each axis 493 LOGICAL , INTENT(in ) , OPTIONAL :: lzfirst ! Whether array being read has been 494 ! ftrans'd to make z index first 478 495 ! 479 496 IF( kiomid > 0 ) THEN 480 497 IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom , cdvar , pv_r1d=pvar, & 481 & ktime=ktime, kstart=kstart, kcount=kcount )498 & ktime=ktime, kstart=kstart, kcount=kcount, lzfirst=lzfirst ) 482 499 ENDIF 483 500 END SUBROUTINE iom_g1d 484 501 485 SUBROUTINE iom_g2d( kiomid, kdom, cdvar, pvar, ktime, kstart, kcount )502 SUBROUTINE iom_g2d( kiomid, kdom, cdvar, pvar, ktime, kstart, kcount, lzfirst ) 486 503 INTEGER , INTENT(in ) :: kiomid ! Identifier of the file 487 504 INTEGER , INTENT(in ) :: kdom ! Type of domain to be read … … 491 508 INTEGER , INTENT(in ), DIMENSION(2) , OPTIONAL :: kstart ! start axis position of the reading 492 509 INTEGER , INTENT(in ), DIMENSION(2) , OPTIONAL :: kcount ! number of points in each axis 510 LOGICAL , INTENT(in ) , OPTIONAL :: lzfirst ! Whether array being read has been 511 ! ftrans'd to make z index first 493 512 ! 494 513 IF( kiomid > 0 ) THEN 495 514 IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom , cdvar , pv_r2d=pvar, & 496 & ktime=ktime, kstart=kstart, kcount=kcount )515 & ktime=ktime, kstart=kstart, kcount=kcount, lzfirst=lzfirst ) 497 516 ENDIF 498 517 END SUBROUTINE iom_g2d 499 518 500 SUBROUTINE iom_g3d( kiomid, kdom, cdvar, pvar, ktime, kstart, kcount )519 SUBROUTINE iom_g3d( kiomid, kdom, cdvar, pvar, ktime, kstart, kcount, lzfirst ) 501 520 INTEGER , INTENT(in ) :: kiomid ! Identifier of the file 502 521 INTEGER , INTENT(in ) :: kdom ! Type of domain to be read … … 506 525 INTEGER , INTENT(in ), DIMENSION(3) , OPTIONAL :: kstart ! start axis position of the reading 507 526 INTEGER , INTENT(in ), DIMENSION(3) , OPTIONAL :: kcount ! number of points in each axis 527 LOGICAL , INTENT(in ), OPTIONAL :: lzfirst ! Whether array being read has been 528 ! ftrans'd to make z index first 508 529 ! 509 530 IF( kiomid > 0 ) THEN 510 531 IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom , cdvar , pv_r3d=pvar, & 511 & ktime=ktime, kstart=kstart, kcount=kcount )532 & ktime=ktime, kstart=kstart, kcount=kcount, lzfirst=lzfirst ) 512 533 ENDIF 513 534 END SUBROUTINE iom_g3d … … 516 537 SUBROUTINE iom_get_123d( kiomid, kdom , cdvar , & 517 538 & pv_r1d, pv_r2d, pv_r3d, & 518 & ktime , kstart, kcount )539 & ktime , kstart, kcount, lzfirst ) 519 540 !!----------------------------------------------------------------------- 520 541 !! *** ROUTINE iom_get_123d *** … … 533 554 INTEGER , DIMENSION(:) , INTENT(in ), OPTIONAL :: kstart ! start position of the reading in each axis 534 555 INTEGER , DIMENSION(:) , INTENT(in ), OPTIONAL :: kcount ! number of points to be read in each axis 556 LOGICAL , INTENT(in ), OPTIONAL :: lzfirst ! Whether array being read has been 557 ! ftrans'd to make z index first 535 558 ! 536 559 LOGICAL :: llnoov ! local definition to read overlap … … 556 579 CHARACTER(LEN=1) :: clrankpv, cldmspc ! 557 580 558 #if defined key_z_first581 !#if defined key_z_first 559 582 !! DCSE_NEMO: need a work array to match layout on disk, which is always z-last 560 583 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: wpv_r3d ! copy of pv_r3d with dimensions permuted … … 562 585 INTEGER, DIMENSION(3) :: ishape_pv_r3d ! size of the dimensions of pv_r3d 563 586 INTEGER :: jk ! loop counter 587 !#endif 588 LOGICAL :: lftrans ! DCSE_NEMO: Whether ftrans is 589 ! in effect for the array we're reading 590 !--------------------------------------------------------------------- 591 ! 592 CALL timing_start('iom_get_123d') 593 594 #if defined key_z_first 595 IF( PRESENT(lzfirst) )THEN 596 lftrans = lzfirst 597 ELSE 598 ! DCSE_NEMO: If we're built with key_z_first then we assume array 599 ! being read has been ftrans'd in the calling routine unless told 600 ! otherwise 601 lftrans = .TRUE. 602 END IF 603 #else 604 ! If we've not been built with key_z_first then we effectively ignore 605 ! the lzfirst argument because ftrans can't have been used. 606 lftrans = .FALSE. 564 607 #endif 565 608 566 !---------------------------------------------------------------------567 !568 609 clname = iom_file(kiomid)%name ! esier to read 569 610 clinfo = ' iom_get_123d, file: '//trim(clname)//', var: '//trim(cdvar) … … 650 691 ELSE 651 692 IF( .NOT. PRESENT(pv_r1d) ) THEN ! not a 1D array 652 IF( idom == jpdom_data ) THEN ; istart(1:2) = (/ mig(1), mjg(1) /) ! icnt(1:2) done bel low653 ELSEIF( idom == jpdom_global ) THEN ; istart(1:2) = (/ nimpp , njmpp /) ! icnt(1:2) done bel low693 IF( idom == jpdom_data ) THEN ; istart(1:2) = (/ mig(1), mjg(1) /) ! icnt(1:2) done below 694 ELSEIF( idom == jpdom_global ) THEN ; istart(1:2) = (/ nimpp , njmpp /) ! icnt(1:2) done below 654 695 ENDIF 655 696 ! we do not read the overlap -> we start to read at nldi, nldj … … 683 724 END DO 684 725 685 #if defined key_z_first686 726 !! DCSE_NEMO: Allocate 3d work-array with z-index last 687 !! to match layout on disk 688 IF ( PRESENT(pv_r3d)) THEN727 !! to match layout on disk if ftrans in use 728 IF (lftrans .AND. PRESENT(pv_r3d)) THEN 689 729 ishape_pv_r3d = SHAPE(pv_r3d) 690 730 IF (ishape_pv_r3d(1) /= jpk) THEN … … 693 733 CALL ctl_warn( trim(clinfo), 'beware: possible problem with 3d array, ', ctmp1 ) 694 734 ENDIF 695 ALLOCATE(wpv_r3d(ishape_pv_r3d(2),ishape_pv_r3d(3),ishape_pv_r3d(1)),STAT=istat_wpv_r3d) 735 ! This assumes that the array we're to read has been ftrans'd in the 736 ! calling routine and therefore that it's z/depth index is its 1st 737 ! whereas on disk it will be its 3rd. 738 ALLOCATE(wpv_r3d(ishape_pv_r3d(2),ishape_pv_r3d(3),ishape_pv_r3d(1)),& 739 STAT=istat_wpv_r3d) 696 740 IF (istat_wpv_r3d /= 0) THEN 697 741 CALL ctl_stop( trim(clinfo), 'failed to allocate wpv_r3d' ) 698 742 ENDIF 699 743 ENDIF 700 #endif701 744 702 745 ! check that icnt matches the input array … … 708 751 IF( irankpv == 1 ) ishape(1:1) = SHAPE(pv_r1d) 709 752 IF( irankpv == 2 ) ishape(1:2) = SHAPE(pv_r2d) 710 #if defined key_z_first 711 IF( irankpv == 3 ) ishape(1:3) = SHAPE(wpv_r3d) 712 #else 713 IF( irankpv == 3 ) ishape(1:3) = SHAPE(pv_r3d) 714 #endif 753 IF( irankpv == 3 )THEN 754 IF( lftrans )THEN 755 ishape(1:3) = SHAPE(wpv_r3d) 756 ELSE 757 ishape(1:3) = SHAPE(pv_r3d) 758 END IF 759 END IF 715 760 ctmp1 = 'd' 716 761 ELSE … … 725 770 ! JMM + SM: ugly patch before getting the new version of lib_mpp) 726 771 ! ishape(1:3) = SHAPE(pv_r3d(nldi:nlei,nldj:nlej,:)) ; ctmp1 = 'd(nldi:nlei,nldj:nlej,:)' 727 #if defined key_z_first 728 IF( llnoov ) THEN 729 ishape(1:3)=SHAPE(wpv_r3d(nldi:nlei,nldj:nlej,:)) 730 ctmp1='d(nldi:nlei,nldj:nlej,:)' 772 IF( lftrans )THEN 773 IF( llnoov ) THEN 774 ishape(1:3)=SHAPE(wpv_r3d(nldi:nlei,nldj:nlej,:)) 775 ctmp1='d(nldi:nlei,nldj:nlej,:)' 776 ELSE 777 ishape(1:3)=SHAPE(wpv_r3d(1 :nlci,1 :nlcj,:)) 778 ctmp1='d(1:nlci,1:nlcj,:)' 779 ENDIF 731 780 ELSE 732 ishape(1:3)=SHAPE(wpv_r3d(1 :nlci,1 :nlcj,:)) 733 ctmp1='d(1:nlci,1:nlcj,:)' 734 ENDIF 735 #else 736 IF( llnoov ) THEN 737 ishape(1:3)=SHAPE(pv_r3d(nldi:nlei,nldj:nlej,:)) 738 ctmp1='d(nldi:nlei,nldj:nlej,:)' 739 ELSE 740 ishape(1:3)=SHAPE(pv_r3d(1 :nlci,1 :nlcj,:)) 741 ctmp1='d(1:nlci,1:nlcj,:)' 742 ENDIF 743 #endif 744 ENDIF 745 ENDIF 781 IF( llnoov ) THEN 782 ishape(1:3)=SHAPE(pv_r3d(nldi:nlei,nldj:nlej,:)) 783 ctmp1='d(nldi:nlei,nldj:nlej,:)' 784 ELSE 785 ishape(1:3)=SHAPE(pv_r3d(1 :nlci,1 :nlcj,:)) 786 ctmp1='d(1:nlci,1:nlcj,:)' 787 ENDIF 788 ENDIF ! ftrans in use 789 ENDIF 790 ENDIF ! ipdom == jpdom_unknown 746 791 747 792 DO jl = 1, irankpv … … 771 816 ENDIF 772 817 773 #if defined key_z_first774 818 SELECT CASE (iom_file(kiomid)%iolib) 775 819 CASE (jpioipsl ) 776 IF ( PRESENT(pv_r3d)) THEN820 IF (lftrans) THEN 777 821 CALL iom_ioipsl_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, wpv_r3d ) 778 822 ELSE 779 CALL iom_ioipsl_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d )823 CALL iom_ioipsl_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, pv_r3d ) 780 824 ENDIF 781 825 CASE (jpnf90 ) 782 IF ( PRESENT(pv_r3d)) THEN826 IF (lftrans) THEN 783 827 CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, wpv_r3d ) 784 828 ELSE 785 CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d )829 CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, pv_r3d ) 786 830 ENDIF 787 831 CASE (jprstdimg) 788 IF ( PRESENT(pv_r3d)) THEN832 IF (lftrans) THEN 789 833 CALL iom_rstdimg_get( kiomid, idom, idvar, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, wpv_r3d ) 790 834 ELSE 791 CALL iom_rstdimg_get( kiomid, idom, idvar, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d )835 CALL iom_rstdimg_get( kiomid, idom, idvar, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, pv_r3d ) 792 836 ENDIF 793 837 CASE DEFAULT 794 838 CALL ctl_stop( TRIM(clinfo)//' accepted IO library are only jpioipsl, jpnf90 and jprstdimg' ) 795 839 END SELECT 796 #else 797 SELECT CASE (iom_file(kiomid)%iolib) 798 CASE (jpioipsl ) ; CALL iom_ioipsl_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, & 799 & pv_r1d, pv_r2d, pv_r3d ) 800 CASE (jpnf90 ) ; CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, & 801 & pv_r1d, pv_r2d, pv_r3d ) 802 CASE (jprstdimg) ; CALL iom_rstdimg_get( kiomid, idom, idvar, ix1, ix2, iy1, iy2, & 803 & pv_r1d, pv_r2d, pv_r3d ) 804 CASE DEFAULT 805 CALL ctl_stop( TRIM(clinfo)//' accepted IO library are only jpioipsl, jpnf90 and jprstdimg' ) 806 END SELECT 807 #endif 808 809 #if defined key_z_first 840 810 841 !! DCSE_NEMO: if necessary, copy 3d work array back into pv_r3d, 811 842 !! and de-allocate the work array 812 IF (PRESENT(pv_r3d)) THEN 813 ! This assumes that pv_r3d is not ftransed 814 DO jk = 1, ishape_pv_r3d(3) 815 DO jj = 1, ishape_pv_r3d(2) 816 DO ji = 1, ishape_pv_r3d(1) 843 IF (lftrans .AND. PRESENT(pv_r3d)) THEN 844 ! This assumes that pv_r3d is ftrans'd in the calling routine so that 845 ! its first dimension rather than last dimension is 'z' 846 ! wpv_r3d is allocated with SHAPE(ishape_pv_r3d(2),ishape_pv_r3d(3),ishape_pv_r3d(1)) 847 DO ji = 1, ishape_pv_r3d(2) 848 DO jj = 1, ishape_pv_r3d(3) 849 DO jk = 1, ishape_pv_r3d(1) 817 850 pv_r3d(jk, ji, jj) = wpv_r3d(ji, jj, jk) 818 851 ENDDO … … 821 854 DEALLOCATE(wpv_r3d) 822 855 ENDIF 823 #endif824 856 825 857 IF( istop == nstop ) THEN ! no additional errors until this point... … … 861 893 ENDIF 862 894 ! 895 CALL timing_stop('iom_get_123d') 896 863 897 END SUBROUTINE iom_get_123d 864 898 … … 883 917 !--------------------------------------------------------------------- 884 918 ! 919 CALL timing_start('iom_gettime') 920 885 921 IF ( PRESENT(cdvar) ) THEN 886 922 tname = cdvar … … 924 960 ENDIF 925 961 ! 962 CALL timing_stop('iom_gettime') 963 ! 926 964 END SUBROUTINE iom_gettime 927 965 … … 961 999 INTEGER :: ivid ! variable id 962 1000 IF( kiomid > 0 ) THEN 1001 ! CALL timing_start('iom_rp0d') 963 1002 IF( iom_file(kiomid)%nfid > 0 ) THEN 964 1003 ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. ) … … 971 1010 END SELECT 972 1011 ENDIF 1012 ! CALL timing_stop('iom_rp0d') 973 1013 ENDIF 974 1014 END SUBROUTINE iom_rp0d … … 983 1023 INTEGER :: ivid ! variable id 984 1024 IF( kiomid > 0 ) THEN 1025 ! CALL timing_start('iom_rp1d') 985 1026 IF( iom_file(kiomid)%nfid > 0 ) THEN 986 1027 ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. ) … … 993 1034 END SELECT 994 1035 ENDIF 1036 ! CALL timing_stop('iom_rp1d') 995 1037 ENDIF 996 1038 END SUBROUTINE iom_rp1d … … 1005 1047 INTEGER :: ivid ! variable id 1006 1048 IF( kiomid > 0 ) THEN 1049 ! CALL timing_start('iom_rp2d') 1007 1050 IF( iom_file(kiomid)%nfid > 0 ) THEN 1008 1051 ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. ) … … 1015 1058 END SELECT 1016 1059 ENDIF 1060 ! CALL timing_stop('iom_rp2d') 1017 1061 ENDIF 1018 1062 END SUBROUTINE iom_rp2d … … 1032 1076 INTEGER :: ji, jj, jk ! Dummy loop indices 1033 1077 IF( kiomid > 0 ) THEN 1078 ! CALL timing_start('iom_rp3d') 1034 1079 IF( iom_file(kiomid)%nfid > 0 ) THEN 1035 1080 ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. ) … … 1056 1101 DEALLOCATE( pvar_trans ) 1057 1102 ENDIF 1103 ! CALL timing_stop('iom_rp3d') 1058 1104 ENDIF 1059 1105 #else 1060 1106 IF( kiomid > 0 ) THEN 1107 ! CALL timing_start('iom_rp3d') 1061 1108 IF( iom_file(kiomid)%nfid > 0 ) THEN 1062 1109 ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. ) … … 1069 1116 END SELECT 1070 1117 ENDIF 1118 ! CALL timing_stop('iom_rp3d') 1071 1119 ENDIF 1072 1120 #endif -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90
r3211 r3432 16 16 !!---------------------------------------------------------------------- 17 17 USE lib_mpp ! distributed memory computing library 18 USE exchmod ! Comms for irregular domain decomposition 18 19 19 20 INTERFACE lbc_lnk 21 #if defined key_mpp_rkpart 22 MODULE PROCEDURE lbc_exch2, lbc_exch3 !, lbc_exch2i, lbc_exch3i 23 #else 20 24 MODULE PROCEDURE mpp_lnk_3d_gather, mpp_lnk_3d, mpp_lnk_2d 25 #endif 21 26 END INTERFACE 22 27 -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
r3211 r3432 65 65 PUBLIC mpp_lnk_3d, mpp_lnk_3d_gather, mpp_lnk_2d, mpp_lnk_2d_e 66 66 PUBLIC mppobc, mpp_ini_ice, mpp_ini_znl 67 PUBLIC mppsize 67 PUBLIC mppsize, MAX_FACTORS, nxfactors, xfactors, nyfactors, yfactors 68 68 PUBLIC lib_mpp_alloc ! Called in nemogcm.F90 69 69 … … 146 146 LOGICAL :: l_isend = .FALSE. ! isend use indicator (T if cn_mpi_send='I') 147 147 INTEGER :: nn_buffer = 0 ! size of the buffer in case of mpi_bsend 148 148 CHARACTER(len=256) :: nn_xfactors = '' ! String holding factors to use for PE grid in x direction 149 CHARACTER(len=256) :: nn_yfactors = '' ! String holding factors to use for PE grid in y direction 150 INTEGER, PARAMETER :: MAX_FACTORS = 20 ! Maximum no. of factors factor() can return 151 ! Arrays to hold specific factorisation of the processor grid specified 152 ! in the namelist. Set to -1 if no specific factorisation requested. 153 INTEGER, SAVE, DIMENSION(MAX_FACTORS) :: xfactors, yfactors 154 INTEGER, SAVE :: nxfactors, nyfactors 155 LOGICAL, PUBLIC :: nn_pttrim = .FALSE. ! Whether to trim dry 156 ! land from PE domains 157 INTEGER, SAVE, PUBLIC :: nn_cpnode = 4 ! Number of cores per 158 ! compute node on current computer 159 149 160 REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: tampon ! buffer in case of bsend 150 161 … … 203 214 !!---------------------------------------------------------------------- 204 215 ! 205 ALLOCATE( t4ns(jpi,jprecj,jpk,2,2) , t4sn(jpi,jprecj,jpk,2,2) , & 216 ALLOCATE( & 217 #if !defined key_mpp_rkpart 218 t4ns(jpi,jprecj,jpk,2,2) , t4sn(jpi,jprecj,jpk,2,2) , & 206 219 & t4ew(jpj,jpreci,jpk,2,2) , t4we(jpj,jpreci,jpk,2,2) , & 207 220 & t4p1(jpi,jprecj,jpk,2,2) , t4p2(jpi,jprecj,jpk,2,2) , & … … 212 225 & t2ew(jpj,jpreci ,2) , t2we(jpj,jpreci ,2) , & 213 226 & t2p1(jpi,jprecj ,2) , t2p2(jpi,jprecj ,2) , & 227 #endif 214 228 ! 215 229 & tr2ns(1-jpr2di:jpi+jpr2di,jprecj+jpr2dj,2) , & … … 248 262 LOGICAL :: mpi_was_called 249 263 ! 250 NAMELIST/nammpp/ cn_mpi_send, nn_buffer, jpni, jpnj, jpnij 264 NAMELIST/nammpp/ cn_mpi_send, nn_buffer, jpni, jpnj, jpnij, & 265 nn_xfactors, nn_yfactors, nn_pttrim, nn_cpnode 251 266 !!---------------------------------------------------------------------- 252 267 ! … … 263 278 WRITE(ldtxt(ii),*) ' mpi send type cn_mpi_send = ', cn_mpi_send ; ii = ii + 1 264 279 WRITE(ldtxt(ii),*) ' size in bytes of exported buffer nn_buffer = ', nn_buffer ; ii = ii + 1 265 280 WRITE(ldtxt(ii),*) ' whether to trim dry points nn_pttrim = ', nn_pttrim ; ii = ii + 1 281 WRITE(ldtxt(ii),*) ' number of cores per compute node nn_cpn = ', nn_cpnode ; ii = ii + 1 266 282 #if defined key_agrif 267 283 IF( .NOT. Agrif_Root() ) THEN … … 284 300 WRITE(ldtxt(ii),*) ' processor grid extent in j jpnj = ',jpnj; ii = ii + 1 285 301 WRITE(ldtxt(ii),*) ' number of local domains jpnij = ',jpnij; ii = ii +1 302 END IF 303 304 ! Check to see whether a specific factorisation of the number of 305 ! processors has been specified in the namelist file 306 nxfactors = 0; xfactors(:) = -1 307 nyfactors = 0; yfactors(:) = -1 308 IF((LEN_TRIM(nn_xfactors) > 0) .OR. (LEN_TRIM(nn_yfactors) > 0))THEN 309 310 IF( (VERIFY(TRIM(nn_xfactors),'0123456789,') > 0) .OR. & 311 (VERIFY(TRIM(nn_yfactors),'0123456789,') > 0) )THEN 312 WRITE(ldtxt(ii),*)'Invalid character in nn_xfactors/nn_yfactors namelist string. '; ii = ii + 1 313 WRITE(ldtxt(ii),*)'Will ignore requested factorisation.' ; ii = ii + 1 314 ELSE 315 READ(nn_xfactors, *,end=80,err=100) xfactors 316 80 CONTINUE 317 READ(nn_yfactors, *,end=90,err=100) yfactors 318 90 CONTINUE 319 320 trim_xarray: DO ji=MAX_FACTORS,1,-1 321 IF (xfactors(ji) .GE. 0) THEN 322 nxfactors = ji 323 EXIT trim_xarray 324 ENDIF 325 ENDDO trim_xarray 326 trim_yarray: DO ji=MAX_FACTORS,1,-1 327 IF (yfactors(ji) .GE. 0) THEN 328 nyfactors = ji 329 EXIT trim_yarray 330 ENDIF 331 ENDDO trim_yarray 332 333 100 CONTINUE 334 WRITE (*,*) 'ARPDBG: n{x,y}factors = ',nxfactors,nyfactors 335 IF(nxfactors < 1 .AND. nyfactors < 1)THEN 336 WRITE(ldtxt(ii),*)'Failed to parse factorisation string' ; ii = ii + 1 337 WRITE(ldtxt(ii),*)' - will ignore requested factorisation.' ; ii = ii + 1 338 ELSE 339 WRITE(ldtxt(ii),*)' automatic factorisation overridden' ; ii = ii + 1 340 WRITE(ldtxt(ii),*)' factors:', xfactors(1:nxfactors), & 341 '-',yfactors(1:nyfactors) 342 ii = ii + 1 343 END IF 344 ENDIF 345 286 346 END IF 287 347 … … 356 416 mynode = mpprank 357 417 ! 418 IF( (jpni < 1) .OR. (jpnj < 1) )THEN 419 IF(mynode == 0)WRITE(ldtxt(ii),*) ' Running on ',mppsize,' MPI processes'; ii = ii + 1 420 END IF 421 358 422 #if defined key_mpp_rep 359 423 CALL MPI_OP_CREATE(DDPDD_MPI, .TRUE., MPI_SUMDD, ierr) … … 363 427 364 428 365 SUBROUTINE mpp_lnk_3d( ptab3d, cd_type, psgn, cd_mpp, pval )429 SUBROUTINE mpp_lnk_3d( ptab3d, cd_type, psgn, cd_mpp, pval, lzero ) 366 430 !!---------------------------------------------------------------------- 367 431 !! *** routine mpp_lnk_3d *** … … 394 458 CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only 395 459 REAL(wp) , OPTIONAL , INTENT(in ) :: pval ! background value (used at closed boundaries) 460 LOGICAL , OPTIONAL , INTENT(in ) :: lzero ! Whether to zero field at closed boundaries 396 461 !! 397 462 INTEGER :: ji, jj, jk, jl ! dummy loop indices … … 400 465 REAL(wp) :: zland 401 466 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat ! for key_mpi_isend 402 !!---------------------------------------------------------------------- 467 LOGICAL :: lzeroarg 468 !!---------------------------------------------------------------------- 469 470 #if defined key_mpp_rkpart 471 CALL ctl_stop('mpp_lnk_3d: should not have been called when key_mpp_rkpart defined!') 472 RETURN 473 #endif 474 ! Deal with optional routine arguments 475 lzeroarg = .TRUE. 476 IF( PRESENT(lzero) ) lzeroarg = lzero 403 477 404 478 IF( PRESENT( pval ) ) THEN ; zland = pval ! set land value … … 443 517 ptab3d(jpi,:,:) = ptab3d( 2 ,:,:) 444 518 ELSE !* closed 445 IF( .NOT. cd_type == 'F' ) ptab3d( 1 :jpreci,:,:) = zland ! south except F-point 446 ptab3d(nlci-jpreci+1:jpi ,:,:) = zland ! north 519 IF( lzeroarg )THEN 520 IF( .NOT. cd_type == 'F' ) ptab3d( 1 :jpreci,:,:) = zland ! south except F-point 521 ptab3d(nlci-jpreci+1:jpi ,:,:) = zland ! north 522 END IF 447 523 ENDIF 448 524 ! ! North-South boundaries (always closed) 449 IF( .NOT. cd_type == 'F' ) ptab3d(:, 1 :jprecj,:) = zland ! south except F-point 450 ptab3d(:,nlcj-jprecj+1:jpj ,:) = zland ! north 525 IF( lzeroarg )THEN 526 IF( .NOT. cd_type == 'F' ) ptab3d(:, 1 :jprecj,:) = zland ! south except F-point 527 ptab3d(:,nlcj-jprecj+1:jpj ,:) = zland ! north 528 END IF 451 529 ! 452 530 ENDIF … … 574 652 575 653 576 SUBROUTINE mpp_lnk_2d( pt2d, cd_type, psgn, cd_mpp, pval )654 SUBROUTINE mpp_lnk_2d( pt2d, cd_type, psgn, cd_mpp, pval, lzero ) 577 655 !!---------------------------------------------------------------------- 578 656 !! *** routine mpp_lnk_2d *** … … 600 678 CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only 601 679 REAL(wp) , OPTIONAL , INTENT(in ) :: pval ! background value (used at closed boundaries) 680 LOGICAL , OPTIONAL , INTENT(in ) :: lzero ! Whether to zero field at closed boundaries 602 681 !! 603 682 INTEGER :: ji, jj, jl ! dummy loop indices … … 606 685 REAL(wp) :: zland 607 686 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat ! for key_mpi_isend 608 !!---------------------------------------------------------------------- 687 LOGICAL :: lzeroarg 688 !!---------------------------------------------------------------------- 689 690 #if defined key_mpp_rkpart 691 CALL ctl_stop('mpp_lnk_3d: should not have been called when key_mpp_rkpart defined!') 692 RETURN 693 #endif 694 695 ! Deal with optional routine arguments 696 lzeroarg = .TRUE. 697 IF( PRESENT(lzero) ) lzeroarg = lzero 609 698 610 699 IF( PRESENT( pval ) ) THEN ; zland = pval ! set land value … … 637 726 pt2d(jpi,:) = pt2d( 2 ,:) ! east 638 727 ELSE ! closed 639 IF( .NOT. cd_type == 'F' ) pt2d( 1 :jpreci,:) = zland ! south except F-point 640 pt2d(nlci-jpreci+1:jpi ,:) = zland ! north 728 IF( lzeroarg )THEN 729 IF( .NOT. cd_type == 'F' ) pt2d( 1 :jpreci,:) = zland ! south except F-point 730 pt2d(nlci-jpreci+1:jpi ,:) = zland ! north 731 END IF 641 732 ENDIF 642 733 ! ! North-South boundaries (always closed) 643 IF( .NOT. cd_type == 'F' ) pt2d(:, 1 :jprecj) = zland !south except F-point 644 pt2d(:,nlcj-jprecj+1:jpj ) = zland ! north 734 IF( lzeroarg )THEN 735 IF( .NOT. cd_type == 'F' ) pt2d(:, 1 :jprecj) = zland !south except F-point 736 pt2d(:,nlcj-jprecj+1:jpj ) = zland ! north 737 END IF 645 738 ! 646 739 ENDIF … … 1782 1875 CALL mppsync 1783 1876 CALL mpi_finalize( info ) 1877 STOP 1784 1878 ! 1785 1879 END SUBROUTINE mppstop -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90
r3211 r3432 129 129 !! 8.5 ! 02-08 (G. Madec) F90 : free form 130 130 !!---------------------------------------------------------------------- 131 USE exchtestmod, ONLY: mpp_test_comms 131 132 INTEGER :: ji, jj, jn ! dummy loop indices 132 133 INTEGER :: ii, ij, ifreq, il1, il2 ! local integers … … 134 135 REAL(wp) :: zidom, zjdom ! local scalars 135 136 INTEGER, DIMENSION(jpni,jpnj) :: iimppt, ijmppt, ilcit, ilcjt ! local workspace 137 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: imask ! Local fake global land mask 136 138 !!---------------------------------------------------------------------- 137 139 … … 424 426 CALL mpp_init_ioipsl 425 427 428 ! ARPDBG - test comms setup 429 ALLOCATE(imask(jpiglo,jpjglo)) 430 imask(:,:) = 1 431 CALL mpp_test_comms(imask) 432 DEALLOCATE(imask) 433 426 434 END SUBROUTINE mpp_init 427 435 -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r3258 r3432 21 21 !! ldf_slp_init : initialization of the slopes computation 22 22 !!---------------------------------------------------------------------- 23 USE timing, ONLY: timing_start, timing_stop 23 24 USE oce ! ocean dynamics and tracers 24 25 USE dom_oce ! ocean space and time domain … … 153 154 REAL(wp) :: zck, zfk, zbw ! - - 154 155 !!---------------------------------------------------------------------- 156 157 CALL timing_start('ldf_slp') 155 158 156 159 IF( wrk_in_use(3, 1) ) THEN … … 475 478 IF( wrk_not_released(3, 1) ) CALL ctl_stop('ldf_slp: failed to release workspace arrays') 476 479 ! 480 CALL timing_stop('ldf_slp', 'section') 481 477 482 END SUBROUTINE ldf_slp 478 483 … … 856 861 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 857 862 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 858 zbu = MIN( zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,ik )* ABS( zau ) )859 zbv = MIN( zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ik )* ABS( zav ) )863 zbu = MIN( zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) ) 864 zbv = MIN( zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) ) 860 865 ! !- Slope at u- & v-points (uslpml, vslpml) 861 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,ik )862 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ik )866 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 867 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 863 868 ! 864 869 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90
r3211 r3432 1251 1251 WRITE(*,*) 'obc_dta: You should not have seen this print! error?', kt 1252 1252 END SUBROUTINE obc_dta 1253 1254 1255 SUBROUTINE obc_dta_bt( kt, kbt ) ! Dummy routine 1256 INTEGER, INTENT( in ) :: kt ! ocean time-step index 1257 INTEGER, INTENT( in ) :: kbt ! barotropic ocean time-step index 1258 WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', kt, kbt 1259 END SUBROUTINE obc_dta_bt 1253 1260 #endif 1254 1261 !!============================================================================== -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90
r3211 r3432 332 332 !!--------------------------------------------------------------------------------- 333 333 CONTAINS 334 SUBROUTINE obc_vol ! Empty routine 334 SUBROUTINE obc_vol(kt) ! Empty routine 335 !! * Arguments 336 INTEGER, INTENT( in ) :: kt 335 337 END SUBROUTINE obc_vol 336 338 #endif -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r3211 r3432 826 826 ! 827 827 IF( wrk_in_use(2, 1) .OR. iwrk_in_use(2,1) ) THEN 828 CALL ctl_stop('fld_weight s: requested workspace arrays are unavailable') ; RETURN828 CALL ctl_stop('fld_weight: requested workspace arrays are unavailable') ; RETURN 829 829 ENDIF 830 830 ! 831 831 IF( nxt_wgt > tot_wgts ) THEN 832 CALL ctl_stop("fld_weight s: weights array size exceeded, increase tot_wgts")832 CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts") 833 833 ENDIF 834 834 ! … … 940 940 941 941 IF( wrk_not_released(2, 1) .OR. & 942 iwrk_not_released(2, 1) ) CALL ctl_stop('fld_weight s: failed to release workspace arrays')942 iwrk_not_released(2, 1) ) CALL ctl_stop('fld_weight: failed to release workspace arrays') 943 943 ! 944 944 END SUBROUTINE fld_weight -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
r2715 r3432 57 57 58 58 # if defined key_lim3 59 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tatm_ice!: air temperature59 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tatm_ice !: air temperature 60 60 # endif 61 61 … … 77 77 & utau_ice(jpi,jpj) , vtau_ice(jpi,jpj) , & 78 78 & fr1_i0 (jpi,jpj) , fr2_i0 (jpi,jpj) , & 79 # if defined key_lim3 80 & emp_ice(jpi,jpj) , tatm_ice(jpi,jpj) , STAT= sbc_ice_alloc ) 81 # else 79 82 & emp_ice(jpi,jpj) , STAT= sbc_ice_alloc ) 83 # endif 80 84 ! 81 85 IF( lk_mpp ) CALL mpp_sum ( sbc_ice_alloc ) -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
r3211 r3432 182 182 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 183 183 ! 184 #if defined key_lim3185 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) !RB ugly patch186 #endif187 ! ! surface ocean fluxes computed with CLIO bulk formulea188 184 IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce_clio( sf, sst_m ) 189 185 ! … … 480 476 zpatm = 101000. ! atmospheric pressure (assumed constant here) 481 477 478 #if defined key_lim3 479 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) ! LIM3: make Tair available in sea-ice. WARNING allocated after call to ice_init 480 #endif 481 ! ! surface ocean fluxes computed with CLIO bulk formulea 482 482 !------------------------------------! 483 483 ! momentum fluxes (utau, vtau ) ! -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r3211 r3432 185 185 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 186 186 187 #if defined key_lim3188 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) ! LIM3: make Tair available in sea-ice189 #endif190 187 ! ! surface ocean fluxes computed with CLIO bulk formulea 191 188 IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce_core( sf, sst_m, ssu_m, ssv_m ) … … 495 492 !!gm end 496 493 494 #if defined key_lim3 495 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) ! LIM3: make Tair available in sea-ice. WARNING allocated after call to ice_init 496 #endif 497 497 ! ----------------------------------------------------------------------------- ! 498 498 ! Wind components and module relative to the moving ocean ( U10m - U_ice ) ! … … 790 790 !!---------------------------------------------------------------------- 791 791 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released, iwrk_in_use, iwrk_not_released 792 USE wrk_nemo, ONLY: dU10 => wrk_2d_1 ! dU [m/s]793 USE wrk_nemo, ONLY: dT => wrk_2d_ 2! air/sea temperature difference [K]794 USE wrk_nemo, ONLY: dq => wrk_2d_ 3! air/sea humidity difference [K]795 USE wrk_nemo, ONLY: Cd_n10 => wrk_2d_ 4! 10m neutral drag coefficient796 USE wrk_nemo, ONLY: Ce_n10 => wrk_2d_ 5! 10m neutral latent coefficient797 USE wrk_nemo, ONLY: Ch_n10 => wrk_2d_ 6! 10m neutral sensible coefficient798 USE wrk_nemo, ONLY: sqrt_Cd_n10 => wrk_2d_ 7! root square of Cd_n10799 USE wrk_nemo, ONLY: sqrt_Cd => wrk_2d_ 8! root square of Cd800 USE wrk_nemo, ONLY: T_vpot => wrk_2d_ 9! virtual potential temperature [K]801 USE wrk_nemo, ONLY: T_star => wrk_2d_ 10! turbulent scale of tem. fluct.802 USE wrk_nemo, ONLY: q_star => wrk_2d_ 11! turbulent humidity of temp. fluct.803 USE wrk_nemo, ONLY: U_star => wrk_2d_ 12! turb. scale of velocity fluct.804 USE wrk_nemo, ONLY: L => wrk_2d_ 13! Monin-Obukov length [m]805 USE wrk_nemo, ONLY: zeta_u => wrk_2d_ 14! stability parameter at height zu806 USE wrk_nemo, ONLY: zeta_t => wrk_2d_ 15! stability parameter at height zt807 USE wrk_nemo, ONLY: U_n10 => wrk_2d_ 16! neutral wind velocity at 10m [m]808 USE wrk_nemo, ONLY: xlogt => wrk_2d_ 17, xct => wrk_2d_18, zpsi_hu => wrk_2d_19, zpsi_ht => wrk_2d_20, zpsi_m => wrk_2d_21792 USE wrk_nemo, ONLY: dU10 => wrk_2d_14 ! dU [m/s] 793 USE wrk_nemo, ONLY: dT => wrk_2d_15 ! air/sea temperature difference [K] 794 USE wrk_nemo, ONLY: dq => wrk_2d_16 ! air/sea humidity difference [K] 795 USE wrk_nemo, ONLY: Cd_n10 => wrk_2d_17 ! 10m neutral drag coefficient 796 USE wrk_nemo, ONLY: Ce_n10 => wrk_2d_18 ! 10m neutral latent coefficient 797 USE wrk_nemo, ONLY: Ch_n10 => wrk_2d_19 ! 10m neutral sensible coefficient 798 USE wrk_nemo, ONLY: sqrt_Cd_n10 => wrk_2d_20 ! root square of Cd_n10 799 USE wrk_nemo, ONLY: sqrt_Cd => wrk_2d_21 ! root square of Cd 800 USE wrk_nemo, ONLY: T_vpot => wrk_2d_22 ! virtual potential temperature [K] 801 USE wrk_nemo, ONLY: T_star => wrk_2d_23 ! turbulent scale of tem. fluct. 802 USE wrk_nemo, ONLY: q_star => wrk_2d_24 ! turbulent humidity of temp. fluct. 803 USE wrk_nemo, ONLY: U_star => wrk_2d_25 ! turb. scale of velocity fluct. 804 USE wrk_nemo, ONLY: L => wrk_2d_26 ! Monin-Obukov length [m] 805 USE wrk_nemo, ONLY: zeta_u => wrk_2d_27 ! stability parameter at height zu 806 USE wrk_nemo, ONLY: zeta_t => wrk_2d_28 ! stability parameter at height zt 807 USE wrk_nemo, ONLY: U_n10 => wrk_2d_29 ! neutral wind velocity at 10m [m] 808 USE wrk_nemo, ONLY: xlogt => wrk_2d_30, xct => wrk_2d_31, zpsi_hu => wrk_2d_32, zpsi_ht => wrk_2d_33, zpsi_m => wrk_2d_34 809 809 USE wrk_nemo, ONLY: stab => iwrk_2d_1 ! 1st guess stability test integer 810 810 !! … … 833 833 !! * Start 834 834 835 IF( wrk_in_use(2, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21) .OR. & 836 iwrk_in_use(2, 1) ) THEN 835 IF( wrk_in_use(2, 14,15,16,17,18,19, & 836 20,21,22,23,24,25,26,27,28,29, & 837 30,31,32,33,34) .OR. & 838 iwrk_in_use(2, 1) ) THEN 837 839 CALL ctl_stop('TURB_CORE_2Z: requested workspace arrays unavailable') ; RETURN 838 END 840 ENDIF 839 841 840 842 !! Initial air/sea differences … … 911 913 END DO 912 914 !! 913 IF( wrk_not_released(2, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21) .OR. & 914 iwrk_not_released(2, 1) ) CALL ctl_stop('TURB_CORE_2Z: requested workspace arrays unavailable') 915 IF( wrk_not_released(2, 14,15,16,17,18,19, & 916 & 20,21,22,23,24,25,26,27,28,29, & 917 & 30,31,32,33,34 ) .OR. & 918 iwrk_not_released(2, 1) ) & 919 CALL ctl_stop('TURB_CORE_2Z: failed to release workspace arrays') 915 920 ! 916 921 END SUBROUTINE TURB_CORE_2Z … … 920 925 !------------------------------------------------------------------------------- 921 926 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 922 USE wrk_nemo, ONLY: X2 => wrk_2d_3 3923 USE wrk_nemo, ONLY: X => wrk_2d_3 4924 USE wrk_nemo, ONLY: stabit => wrk_2d_3 5927 USE wrk_nemo, ONLY: X2 => wrk_2d_35 928 USE wrk_nemo, ONLY: X => wrk_2d_36 929 USE wrk_nemo, ONLY: stabit => wrk_2d_37 925 930 !! 926 931 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta … … 930 935 !------------------------------------------------------------------------------- 931 936 932 IF( wrk_in_use(2, 3 3,34,35) ) THEN937 IF( wrk_in_use(2, 35,36,37) ) THEN 933 938 CALL ctl_stop('psi_m: requested workspace arrays unavailable') ; RETURN 934 939 ENDIF … … 939 944 & + (1. - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable 940 945 941 IF( wrk_not_released(2, 3 3,34,35) ) CALL ctl_stop('psi_m: failed to release workspace arrays')946 IF( wrk_not_released(2, 35,36,37) ) CALL ctl_stop('psi_m: failed to release workspace arrays') 942 947 ! 943 948 END FUNCTION psi_m … … 947 952 !------------------------------------------------------------------------------- 948 953 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 949 USE wrk_nemo, ONLY: X2 => wrk_2d_3 3950 USE wrk_nemo, ONLY: X => wrk_2d_3 4951 USE wrk_nemo, ONLY: stabit => wrk_2d_3 5954 USE wrk_nemo, ONLY: X2 => wrk_2d_35 955 USE wrk_nemo, ONLY: X => wrk_2d_36 956 USE wrk_nemo, ONLY: stabit => wrk_2d_37 952 957 ! 953 958 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta … … 956 961 !------------------------------------------------------------------------------- 957 962 958 IF( wrk_in_use(2, 3 3,34,35) ) THEN963 IF( wrk_in_use(2, 35,36,37) ) THEN 959 964 CALL ctl_stop('psi_h: requested workspace arrays unavailable') ; RETURN 960 965 ENDIF … … 965 970 & + (1. - stabit)*(2.*log( (1. + X2)/2. )) ! Unstable 966 971 967 IF( wrk_not_released(2, 3 3,34,35) ) CALL ctl_stop('psi_h: failed to release workspace arrays')972 IF( wrk_not_released(2, 35,36,37) ) CALL ctl_stop('psi_h: failed to release workspace arrays') 968 973 ! 969 974 END FUNCTION psi_h -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r3211 r3432 95 95 !!--------------------------------------------------------------------- 96 96 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 97 USE wrk_nemo, ONLY: zalb_ice_os => wrk_3d_1 ! albedo of the ice under overcast sky 98 USE wrk_nemo, ONLY: zalb_ice_cs => wrk_3d_2 ! albedo of ice under clear sky 97 USE wrk_nemo, ONLY: wrk_3d_1, wrk_3d_2 ! for albedo of ice under overcast/clear sky 99 98 !! 100 99 INTEGER, INTENT(in) :: kt ! ocean time step … … 103 102 INTEGER :: jl ! dummy loop index 104 103 REAL(wp) :: zcoef ! local scalar 104 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_ice_os, zalb_ice_cs ! albedo of the ice under overcast/clear sky 105 105 !!---------------------------------------------------------------------- 106 106 … … 108 108 CALL ctl_stop( 'sbc_ice_lim: requested workspace arrays are unavailable' ) ; RETURN 109 109 ENDIF 110 zalb_ice_os => wrk_3d_1(:,:,1:jpl) ; zalb_ice_cs => wrk_3d_2(:,:,1:jpl) 110 111 111 112 IF( kt == nit000 ) THEN -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
r3211 r3432 50 50 51 51 INTEGER , PUBLIC :: nkrnf = 0 !: nb of levels over which Kz is increased at river mouths 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnfmsk !: river mouth mask (hori.) 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: rnfmsk_z !: river mouth mask (vert.) 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_rnf !: depth of runoff in m 55 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nk_rnf !: depth of runoff in model levels 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rnf_tsc_b, rnf_tsc !: before and now T & S runoff contents 57 ! ! [K.m/s & PSU.m/s] 58 52 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnfmsk !: river mouth mask (hori.) 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: rnfmsk_z !: river mouth mask (vert.) 54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_rnf !: depth of runoff in m 55 INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nk_rnf !: depth of runoff in model levels 56 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rnf_tsc_b, rnf_tsc !: before & now T & S runoff contents [K.m/s & PSU.m/s] 57 ! ARP put PUBLIC qualifier here to bring length of above line below 132 chars 58 PUBLIC :: rnfmsk, rnfmsk_z, h_rnf, nk_rnf, rnf_tsc_b, rnf_tsc 59 59 60 REAL(wp) :: r1_rau0 ! = 1 / rau0 60 61 61 62 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnf ! structure: river runoff (file information, fields read) 62 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_s_rnf ! structure: river runoff salinity (file info , fields read)63 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_t_rnf ! structure: river runoff temperature (file info , fields read)63 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_s_rnf ! structure: river runoff salinity (file information, fields read) 64 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_t_rnf ! structure: river runoff temperature (file information, fields read) 64 65 65 66 !! * Control permutation of array indices … … 447 448 !!---------------------------------------------------------------------- 448 449 ! 449 INTEGER :: inum ! temporary integers450 CHARACTER(len= 32) :: cl_rnfile ! runoff file name450 INTEGER :: inum ! temporary integers 451 CHARACTER(len=140) :: cl_rnfile ! runoff file name 451 452 !!---------------------------------------------------------------------- 452 453 ! -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r3211 r3432 71 71 !! - save the trends 72 72 !!---------------------------------------------------------------------- 73 USE timing, ONLY: timing_start, timing_stop 73 74 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 74 75 USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace … … 108 109 !!---------------------------------------------------------------------- 109 110 111 CALL timing_start('tra_adv_tvd') 112 110 113 IF( wrk_in_use(3, 12,13) ) THEN 111 114 CALL ctl_stop('tra_adv_tvd: requested workspace arrays unavailable') ; RETURN … … 140 143 ! -------------------------------------------------------------------- 141 144 ! upstream tracer flux in the i and j direction 145 CALL timing_start('tvd_upstream') 142 146 #if defined key_z_first 143 147 DO jj = 1, jpjm1 … … 159 163 END DO 160 164 END DO 165 CALL timing_stop('tvd_upstream','section') 161 166 162 167 ! upstream tracer flux in the k direction 163 168 ! Surface value 169 CALL timing_start('tvd_upstreamk') 164 170 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! volume variable 165 171 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface … … 181 187 END DO 182 188 END DO 189 CALL timing_stop('tvd_upstreamk','section') 183 190 184 191 ! total advective trend 192 CALL timing_start('tvd_tot') 185 193 #if defined key_z_first 186 194 DO jj = 2, jpjm1 … … 205 213 END DO 206 214 END DO 215 CALL timing_stop('tvd_tot','section') 216 207 217 ! ! Lateral boundary conditions on zwi (unchanged sign) 218 CALL timing_start('tvd_lbc') 208 219 CALL lbc_lnk( zwi, 'T', 1. ) 220 CALL timing_stop('tvd_lbc','section') 209 221 210 222 ! ! trend diagnostics (contribution of upstream fluxes) … … 222 234 ! -------------------------------------------------- 223 235 ! antidiffusive flux on i and j 236 CALL timing_start('tvd_antidiff') 224 237 #if defined key_z_first 225 238 DO jj = 1, jpjm1 … … 236 249 END DO 237 250 END DO 251 CALL timing_stop('tvd_antidiff','section') 238 252 239 253 ! antidiffusive flux on k 254 CALL timing_start('tvd_antidiffk') 240 255 #if defined key_z_first 241 256 DO jj = 1, jpj … … 254 269 END DO 255 270 END DO 271 CALL timing_stop('tvd_antidiffk','section') 272 273 CALL timing_start('tvd_lbc') 256 274 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions 257 275 CALL lbc_lnk( zwz, 'W', 1. ) 276 CALL timing_stop('tvd_lbc','section') 258 277 259 278 ! 4. monotonicity algorithm 260 279 ! ------------------------- 280 CALL timing_start('tvd_nonosc') 261 281 CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 282 CALL timing_stop('tvd_nonosc','section') 262 283 263 284 264 285 ! 5. final trend with corrected fluxes 265 286 ! ------------------------------------ 287 CALL timing_start('tvd_finaltr') 266 288 #if defined key_z_first 267 289 DO jj = 2, jpjm1 … … 300 322 ENDIF 301 323 ! 324 CALL timing_stop('tvd_finaltr','section') 325 302 326 END DO 303 327 ! … … 308 332 IF( wrk_not_released(3, 12,13) ) CALL ctl_stop('tra_adv_tvd: failed to release workspace arrays') 309 333 ! 334 CALL timing_stop('tra_adv_tvd','section') 310 335 311 336 !! * Reset control of array index permutation … … 334 359 USE wrk_nemo, ONLY: zbetup => wrk_3d_8 , zbetdo => wrk_3d_9 ! 3D workspace 335 360 USE wrk_nemo, ONLY: zbup => wrk_3d_10 , zbdo => wrk_3d_11 ! - - 361 USE timing, ONLY: timing_start, timing_stop 336 362 337 363 !! DCSE_NEMO: need additional directives for renamed module variables … … 418 444 END DO 419 445 END DO 446 CALL timing_start('tvd_lbc') 420 447 CALL lbc_lnk( zbetup, 'T', 1. ) ; CALL lbc_lnk( zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) 448 CALL timing_stop('tvd_lbc','section') 421 449 422 450 … … 452 480 END DO 453 481 END DO 482 CALL timing_start('tvd_lbc') 454 483 CALL lbc_lnk( paa, 'U', -1. ) ; CALL lbc_lnk( pbb, 'V', -1. ) ! lateral boundary condition (changed sign) 484 CALL timing_stop('tvd_lbc','section') 485 455 486 ! 456 487 IF( wrk_not_released(3, 8,9,10,11) ) CALL ctl_stop('nonosc: failed to release workspace arrays') -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
r3211 r3432 98 98 !! ** Action : Update pta arrays with the before rotated diffusion 99 99 !!---------------------------------------------------------------------- 100 USE timing, ONLY: timing_start, timing_stop 100 101 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 101 102 USE oce , ONLY: zftu => ua , zftv => va ! (ua,va) used as workspace … … 103 104 !FTRANS zftu zftv :I :I :z 104 105 #if defined key_z_first 105 USE wrk_nemo, ONLY: wdkt => wrk_3d_9 , wdk1t => wrk_3d_10 ! 3D workspace106 ! USE wrk_nemo, ONLY: wdkt => wrk_3d_9 , wdk1t => wrk_3d_10 ! 3D workspace 106 107 !FTRANS wdkt wdk1t :I :I :z 107 108 #else … … 131 132 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 132 133 REAL(wp) :: zcoef0, zbtr, ztra ! - - 134 #if defined key_z_first 135 REAL(wp) :: wdkt , wdki1t , wdkim1t , wdkj1t , wdkjm1t 136 REAL(wp) :: wdk1t, wdk1i1t, wdk1im1t, wdk1j1t, wdk1jm1t 137 #endif 138 133 139 #if defined key_diaar5 134 140 REAL(wp) :: zztmp ! local scalar 135 141 #endif 136 142 !!---------------------------------------------------------------------- 143 144 CALL timing_start('tra_ldf_iso') 137 145 138 146 #if defined key_z_first … … 151 159 ! 152 160 ! ! =========== 161 !DIR$ SHORTLOOP 153 162 DO jn = 1, kjpt ! tracer loop 154 163 ! ! =========== … … 157 166 !! I - masked horizontal derivative 158 167 !!---------------------------------------------------------------------- 168 CALL timing_start('traldf_iso_I') 159 169 !!bug ajout.... why? ( 1,jpj,:) and (jpi,1,:) should be sufficient.... 170 #if defined key_z_first 171 DO jj=1,jpj,1 172 DO jk=1,jpk,1 173 zdit(1 ,jj,jk) = 0.0_wp 174 zdit(jpi,jj,jk) = 0.0_wp 175 zdjt(1 ,jj,jk) = 0.0_wp 176 zdjt(jpi,jj,jk) = 0.0_wp 177 END DO 178 END DO 179 #else 160 180 zdit (1,:,:) = 0.e0 ; zdit (jpi,:,:) = 0.e0 161 181 zdjt (1,:,:) = 0.e0 ; zdjt (jpi,:,:) = 0.e0 182 #endif 162 183 !!end 163 184 … … 185 206 END DO 186 207 ENDIF 208 ! 209 CALL timing_stop('traldf_iso_I','section') 187 210 188 211 !!---------------------------------------------------------------------- 189 212 !! II - horizontal trend (full) 190 213 !!---------------------------------------------------------------------- 214 CALL timing_start('traldf_iso_II') 191 215 #if defined key_z_first 192 216 ! 1. Vertical tracer gradient at level jk and jk+1 … … 194 218 ! surface boundary condition: wdkt(jk=1)=wdkt(jk=2) 195 219 196 DO jj = 1, jpj197 DO ji = 1, jpi198 DO jk = 1, jpkm1199 wdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1)200 END DO201 wdkt(ji,jj,1) = wdk1t(ji,jj,1)202 DO jk = 2, jpkm1203 wdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk)204 END DO205 END DO206 END DO220 !!$ DO jj = 1, jpj 221 !!$ DO ji = 1, jpi 222 !!$ DO jk = 1, jpkm1 223 !!$ wdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 224 !!$ END DO 225 !!$ wdkt(ji,jj,1) = wdk1t(ji,jj,1) 226 !!$ DO jk = 2, jpkm1 227 !!$ wdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 228 !!$ END DO 229 !!$ END DO 230 !!$ END DO 207 231 208 232 ! 2. Horizontal fluxes 209 233 ! -------------------- 210 DO jj = 1 , jpjm1 211 DO ji = 1, jpim1 212 DO jk = 1, jpkm1 213 zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) 214 zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) 215 zmsku = 1. / MAX( tmask(ji+1,jj,jk ) + tmask(ji,jj,jk+1) & 216 & + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk ), 1. ) 217 zmskv = 1. / MAX( tmask(ji,jj+1,jk ) + tmask(ji,jj,jk+1) & 218 & + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk ), 1. ) 219 zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 220 zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 221 zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & 222 & + zcof1 * ( wdkt (ji+1,jj,jk) + wdk1t(ji,jj,jk) & 223 & + wdk1t(ji+1,jj,jk) + wdkt (ji,jj,jk) ) ) * umask(ji,jj,jk) 224 zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 225 & + zcof2 * ( wdkt (ji,jj+1,jk) + wdk1t(ji,jj,jk) & 226 & + wdk1t(ji,jj+1,jk) + wdkt (ji,jj,jk) ) ) * vmask(ji,jj,jk) 227 END DO 228 END DO 229 END DO 230 231 ! II.4 Second derivative (divergence) and add to the general trend 232 ! ---------------------------------------------------------------- 234 !!$ DO jj = 1 , jpjm1 235 !!$ DO ji = 1, jpim1 236 !!$ DO jk = 1, jpkm1 237 !!$ zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) 238 !!$ zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) 239 !!$ zmsku = 1. / MAX( tmask(ji+1,jj,jk ) + tmask(ji,jj,jk+1) & 240 !!$ & + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk ), 1. ) 241 !!$ zmskv = 1. / MAX( tmask(ji,jj+1,jk ) + tmask(ji,jj,jk+1) & 242 !!$ & + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk ), 1. ) 243 !!$ zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 244 !!$ zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 245 !!$ zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & 246 !!$ & + zcof1 * ( wdkt (ji+1,jj,jk) + wdk1t(ji,jj,jk) & 247 !!$ & + wdk1t(ji+1,jj,jk) + wdkt (ji,jj,jk) ) ) * umask(ji,jj,jk) 248 !!$ zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 249 !!$ & + zcof2 * ( wdkt (ji,jj+1,jk) + wdk1t(ji,jj,jk) & 250 !!$ & + wdk1t(ji,jj+1,jk) + wdkt (ji,jj,jk) ) ) * vmask(ji,jj,jk) 251 !!$ END DO 252 !!$ END DO 253 !!$ END DO 254 233 255 DO jj = 2 , jpjm1 234 256 DO ji = 2, jpim1 235 257 DO jk = 1, jpkm1 236 zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 237 ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) 258 259 ! 1. Vertical tracer gradient at level jk and jk+1 260 ! ------------------------------------------------ 261 ! surface boundary condition: wdkt(jk=1)=wdkt(jk=2) 262 263 wdk1t = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 264 wdk1i1t = ( ptb(ji+1,jj,jk,jn) - ptb(ji+1,jj,jk+1,jn) ) * tmask(ji+1,jj,jk+1) 265 wdk1im1t = ( ptb(ji-1,jj,jk,jn) - ptb(ji-1,jj,jk+1,jn) ) * tmask(ji-1,jj,jk+1) 266 wdk1j1t = ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj+1,jk+1,jn) ) * tmask(ji,jj+1,jk+1) 267 wdk1jm1t = ( ptb(ji,jj-1,jk,jn) - ptb(ji,jj-1,jk+1,jn) ) * tmask(ji,jj-1,jk+1) 268 269 IF(jk > 1)THEN 270 wdkt = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 271 wdki1t = ( ptb(ji+1,jj,jk-1,jn) - ptb(ji+1,jj,jk,jn) ) * tmask(ji+1,jj,jk) 272 wdkim1t = ( ptb(ji-1,jj,jk-1,jn) - ptb(ji-1,jj,jk,jn) ) * tmask(ji-1,jj,jk) 273 wdkj1t = ( ptb(ji,jj+1,jk-1,jn) - ptb(ji,jj+1,jk,jn) ) * tmask(ji,jj+1,jk) 274 wdkjm1t = ( ptb(ji,jj-1,jk-1,jn) - ptb(ji,jj-1,jk,jn) ) * tmask(ji,jj-1,jk) 275 ELSE 276 wdkt = wdk1t 277 wdki1t = wdk1i1t 278 wdkim1t= wdk1im1t 279 wdkj1t = wdk1j1t 280 wdkjm1t= wdk1jm1t 281 END IF 282 283 ! II.4 Second derivative (divergence) and add to the general trend 284 ! ---------------------------------------------------------------- 285 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 286 287 ztra = zbtr * ( & 288 289 ! zftu(ji,jj,jk) - 290 ( ((fsahtu(ji,jj,jk) + pahtb0) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)) * zdit(ji,jj,jk) & 291 - ( fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) / & 292 MAX( tmask(ji+1,jj,jk ) + tmask(ji,jj,jk+1) & 293 + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk ), 1.) ) * & 294 (wdki1t + wdk1t + wdk1i1t + wdkt) ) * umask(ji,jj,jk) - & 295 296 ! zftu(ji-1,jj,jk) + 297 ( ((fsahtu(ji-1,jj,jk) + pahtb0) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) / e1u(ji-1,jj)) * zdit(ji-1,jj,jk) & 298 - ( fsahtu(ji-1,jj,jk) * e2u(ji-1,jj) * uslp(ji-1,jj,jk) / & 299 MAX( tmask(ji,jj,jk ) + tmask(ji-1,jj,jk+1) & 300 + tmask(ji,jj,jk+1) + tmask(ji-1,jj,jk ), 1.) ) * & 301 (wdkt + wdk1im1t + wdk1t + wdkim1t) ) * umask(ji-1,jj,jk) + & 302 303 304 ! zftv(ji,jj,jk) - 305 ( ((fsahtv(ji,jj,jk) + pahtb0) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)) * zdjt(ji,jj,jk) & 306 & - ( fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) / & 307 MAX( tmask(ji,jj+1,jk ) + tmask(ji,jj,jk+1) & 308 & + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk ), 1. )) * & 309 (wdkj1t + wdk1t + wdk1j1t + wdkt) ) * vmask(ji,jj,jk) - & 310 ! zftv(ji,jj-1,jk) & 311 ( ((fsahtv(ji,jj-1,jk) + pahtb0) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) / e2v(ji,jj-1)) * zdjt(ji,jj-1,jk) & 312 & - ( fsahtv(ji,jj-1,jk) * e1v(ji,jj-1) * vslp(ji,jj-1,jk) / & 313 MAX( tmask(ji,jj,jk ) + tmask(ji,jj-1,jk+1) & 314 & + tmask(ji,jj,jk+1) + tmask(ji,jj-1,jk ), 1. )) * & 315 (wdkt + wdk1jm1t + wdk1t + wdkjm1t) ) * vmask(ji,jj-1,jk) & 316 317 ) 238 318 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 239 319 END DO … … 295 375 ! "Poleward" diffusive heat or salt transports (T-S case only) 296 376 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 297 IF( jn == jp_tem) htr_ldf(:) = ptr_vj( zftv(:,:,:) ) 298 IF( jn == jp_sal) str_ldf(:) = ptr_vj( zftv(:,:,:) ) 377 IF( jn == jp_tem)THEN 378 htr_ldf = ptr_vj( zftv ) 379 END IF 380 IF( jn == jp_sal)THEN 381 str_ldf = ptr_vj( zftv ) 382 END IF 299 383 ENDIF 300 384 … … 338 422 END IF 339 423 #endif 424 CALL timing_stop('traldf_iso_II','section') 340 425 341 426 !!---------------------------------------------------------------------- 342 427 !! III - vertical trend of T & S (extra diagonal terms only) 343 428 !!---------------------------------------------------------------------- 429 CALL timing_start('traldf_iso_III') 344 430 345 431 ! Local constant initialization 346 432 ! ----------------------------- 433 #if defined key_z_first 434 DO jj=1,jpj,1 435 DO jk=1,jpk,1 436 ztfw(1 ,jj,jk) = 0.0_wp 437 ztfw(jpi,jj,jk) = 0.0_wp 438 END DO 439 END DO 440 #else 347 441 ztfw(1,:,:) = 0.e0 ; ztfw(jpi,:,:) = 0.e0 348 442 #endif 349 443 ! Vertical fluxes 350 444 ! --------------- 351 445 352 446 ! Surface and bottom vertical fluxes set to zero 447 #if defined key_z_first 448 DO ji=1,jpi,1 449 DO jj=1,jpj,1 450 ztfw(ji,jj,1 ) = 0.0_wp 451 ztfw(ji,jj,jpk) = 0.0_wp 452 END DO 453 END DO 454 #else 353 455 ztfw(:,:, 1 ) = 0.e0 ; ztfw(:,:,jpk) = 0.e0 354 456 #endif 457 355 458 ! interior (2=<jk=<jpk-1) 356 459 #if defined key_z_first … … 400 503 END DO 401 504 ! 505 506 CALL timing_stop('traldf_iso_III','section') 507 402 508 END DO 403 509 ! … … 409 515 wrk_not_released(2, 1,2,3) ) CALL ctl_stop('tra_ldf_iso: failed to release workspace arrays') 410 516 #endif 517 ! 518 CALL timing_stop('tra_ldf_iso','section') 411 519 ! 412 520 END SUBROUTINE tra_ldf_iso -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
r3211 r3432 85 85 !! ** Action : - pta becomes the after tracer 86 86 !!--------------------------------------------------------------------- 87 USE timing, ONLY: timing_start, timing_stop 87 88 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 88 89 USE oce , ONLY: zwd => ua , zws => va ! (ua,va) used as 3D workspace … … 109 110 REAL(wp) :: zrhs, ze3tb, ze3tn, ze3ta ! local scalars 110 111 !!--------------------------------------------------------------------- 112 113 CALL timing_start('tra_zdf_imp') 111 114 112 115 IF( wrk_in_use(3, 6,7) ) THEN … … 202 205 DO ji = 2, jpim1 203 206 DO jk = 1, jpkm1 204 ze3ta = ( 1. - r_vvl ) + r_vvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point 205 ze3tn = r_vvl + ( 1. - r_vvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point 207 ! after scale factor at T-point 208 ze3ta = ( 1. - r_vvl ) + r_vvl * fse3t_a(ji,jj,jk) 209 ! now scale factor at T-point 210 ze3tn = r_vvl + ( 1. - r_vvl ) * fse3t_n(ji,jj,jk) 206 211 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) 207 212 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) … … 330 335 IF( wrk_not_released(3, 6,7) ) CALL ctl_stop('tra_zdf_imp: failed to release workspace arrays') 331 336 ! 337 CALL timing_stop('tra_zdf_imp','section') 338 ! 332 339 END SUBROUTINE tra_zdf_imp 333 340 -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90
r3211 r3432 263 263 #else 264 264 DO jk = 1, jpkm1 265 !! DCSE_NEMO: This looks plain wrong!266 ! t2(ktrd) = t2(ktrd) + SUM( ptrd3dx(ji,jj,jk) * tn(ji,jj,jk) * e1e2t(:,:) * fse3t(:,:,jk) )267 ! s2(ktrd) = s2(ktrd) + SUM( ptrd3dy(ji,jj,jk) * sn(ji,jj,jk) * e1e2t(:,:) * fse3t(:,:,jk) )268 265 t2(ktrd) = t2(ktrd) + SUM( ptrd3dx(:,:,jk) * tn(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) ) 269 266 s2(ktrd) = s2(ktrd) + SUM( ptrd3dy(:,:,jk) * sn(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) ) -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r3211 r3432 61 61 ! !!! ** Namelist namzdf_gls ** 62 62 LOGICAL :: ln_crban = .FALSE. ! =T use Craig and Banner scheme 63 LOGICAL :: ln_length_lim = .FALSE. ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 64 LOGICAL :: ln_sigpsi = .FALSE. ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing 63 LOGICAL :: ln_length_lim = .FALSE. ! use limit on the dissipation rate under 64 ! stable stratification (Galperin et al. 1988) 65 LOGICAL :: ln_sigpsi = .FALSE. ! Activate Burchard (2003) modification for 66 ! k-eps closure & wave breaking mixing 65 67 INTEGER :: nn_tkebc_surf = 0 ! TKE surface boundary condition (=0/1) 66 68 INTEGER :: nn_tkebc_bot = 0 ! TKE bottom boundary condition (=0/1) … … 72 74 REAL(wp) :: rn_epsmin = 1.e-12_wp ! minimum value of dissipation (m2/s3) 73 75 REAL(wp) :: rn_emin = 1.e-6_wp ! minimum value of TKE (m2/s2) 74 REAL(wp) :: rn_charn = 2.e+5_wp ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 76 REAL(wp) :: rn_charn = 2.e+5_wp ! Charnock constant for surface breaking waves 77 ! mixing : 1400. (standard) or 2.e5 (Stacey value) 75 78 REAL(wp) :: rn_crban = 100._wp ! Craig and Banner constant for surface breaking waves mixing 76 79 … … 167 170 USE wrk_nemo, ONLY: eps => wrk_3d_4 ! dissipation rate 168 171 USE wrk_nemo, ONLY: zwall_psi => wrk_3d_5 ! Wall function use in the wb case (ln_sigpsi.AND.ln_crban=T) 172 USE timing, ONLY: timing_start, timing_stop 169 173 170 174 !! DCSE_NEMO: need additional directives for renamed module variables … … 180 184 REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - - 181 185 !!-------------------------------------------------------------------- 186 187 CALL timing_start('zdf_gls') 182 188 183 189 IF( wrk_in_use(2, 1,2,3) .OR. wrk_in_use(3, 1,2,3,4,5) ) THEN … … 319 325 IF( ln_sigpsi ) THEN 320 326 zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1. 321 zwall_psi(ji,jj,jk) = rsc_psi / ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 327 zwall_psi(ji,jj,jk) = rsc_psi / ( zsigpsi * rsc_psi + & 328 (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 322 329 ELSE 323 330 zwall_psi(ji,jj,jk) = 1._wp … … 1048 1055 DO ji = fs_2, fs_jpim1 ! vector opt. 1049 1056 #endif 1050 DO jk = 1, jpk1051 DO jj = 2, jpjm11052 DO ji = fs_2, fs_jpim1 ! vector opt.1053 1057 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * mxln(ji,jj,jk) 1054 1058 zav = zsqen * avmu(ji,jj,jk) … … 1089 1093 IF( wrk_not_released(2, 1,2,3) .OR. & 1090 1094 wrk_not_released(3, 1,2,3,4,5) ) CALL ctl_stop('zdf_gls: failed to release workspace arrays') 1095 ! 1096 CALL timing_stop('zdf_gls','section') 1091 1097 ! 1092 1098 END SUBROUTINE zdf_gls -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
r3211 r3432 16 16 USE iom ! I/O library 17 17 USE lib_mpp ! MPP library 18 USE trc_oce, ONLY : lk_offline ! offline flag 18 19 19 20 IMPLICIT NONE … … 45 46 !! *** FUNCTION zdf_mxl_alloc *** 46 47 !!---------------------------------------------------------------------- 47 ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 48 ! 49 IF( lk_mpp ) CALL mpp_sum ( zdf_mxl_alloc ) 50 IF( zdf_mxl_alloc /= 0 ) CALL ctl_warn('zdf_mxl_alloc: failed to allocate arrays.') 48 zdf_mxl_alloc = 0 ! set to zero if no array to be allocated 49 IF( .NOT. ALLOCATED( nmln ) ) THEN 50 ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 51 ! 52 IF( lk_mpp ) CALL mpp_sum ( zdf_mxl_alloc ) 53 IF( zdf_mxl_alloc /= 0 ) CALL ctl_warn('zdf_mxl_alloc: failed to allocate arrays.') 54 ! 55 ENDIF 51 56 END FUNCTION zdf_mxl_alloc 52 57 … … 124 129 END DO 125 130 END DO 126 CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth 127 CALL iom_put( "mldkz5" , hmld ) ! turbocline depth 131 IF( .NOT.lk_offline ) THEN ! no need to output in offline mode 132 CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth 133 CALL iom_put( "mldkz5" , hmld ) ! turbocline depth 134 ENDIF 128 135 129 136 IF(ln_ctl) CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 ) -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r3187 r3432 73 73 USE mod_ioclient 74 74 #endif 75 USE partition_mod ! irregular domain partitioning 76 USE timing, ONLY: timing_init, timing_finalize, timing_disable, timing_enable 77 78 !#define ARPDEBUG 75 79 76 80 IMPLICIT NONE … … 125 129 IF(lwp) WRITE(numout,cform_aaa) ! Flag AAAAAAA 126 130 131 CALL timing_enable() 127 132 ! !-----------------------! 128 133 ! !== time stepping ==! … … 171 176 ENDIF 172 177 ! 178 CALL timing_finalize ! Timing report 179 173 180 CALL nemo_closefile 174 181 #if defined key_oasis3 || defined key_oasis4 … … 189 196 INTEGER :: ji ! dummy loop indices 190 197 INTEGER :: ilocal_comm ! local integer 191 CHARACTER(len=80), DIMENSION( 16) :: cltxt198 CHARACTER(len=80), DIMENSION(24) :: cltxt 192 199 !! 193 200 NAMELIST/namctl/ ln_ctl , nn_print, nn_ictls, nn_ictle, & … … 195 202 !!---------------------------------------------------------------------- 196 203 ! 197 cltxt = ''204 cltxt(:) = '' 198 205 ! 199 206 ! ! open Namelist file … … 228 235 lwp = (narea == 1) .OR. ln_ctl ! control of all listing output print 229 236 237 CALL timing_init ! Init timing module 238 CALL timing_disable ! but disable during startup 239 230 240 ! If dimensions of processor grid weren't specified in the namelist file 231 241 ! then we calculate them here now that we have our communicator size 232 242 IF( (jpni < 1) .OR. (jpnj < 1) )THEN 233 243 #if defined key_mpp_mpi 244 #if defined key_mpp_rkpart 245 IF( Agrif_Root() ) CALL nemo_recursive_partition(mppsize) 246 #else 234 247 IF( Agrif_Root() ) CALL nemo_partition(mppsize) 248 #endif 235 249 #else 236 250 jpni = 1 … … 244 258 ! than variables 245 259 IF( Agrif_Root() ) THEN 260 #if ! defined key_mpp_rkpart 246 261 jpi = ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci ! first dim. 247 262 jpj = ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj ! second dim. 248 jpk = jpkdta ! third dim249 263 jpim1 = jpi-1 ! inner domain indices 250 264 jpjm1 = jpj-1 ! " " 251 jpkm1 = jpk-1 ! " "252 265 jpij = jpi*jpj ! jpi x j 266 #endif 267 jpk = jpkdta ! third dim 268 jpkm1 = jpk-1 ! inner domain indices 253 269 ENDIF 254 270 … … 264 280 WRITE(numout,*) 265 281 WRITE(numout,*) 266 DO ji = 1, SIZE(cltxt )282 DO ji = 1, SIZE(cltxt,1) 267 283 IF( TRIM(cltxt(ji)) /= '' ) WRITE(numout,*) cltxt(ji) ! control print of mynode 268 284 END DO … … 282 298 283 299 ! ! Domain decomposition 300 #if defined key_mpp_rkpart 301 CALL mpp_init3 ! Remainder of set-up for 302 ! recursive partitioning 303 #else 284 304 IF( jpni*jpnj == jpnij ) THEN ; CALL mpp_init ! standard cutting out 285 305 ELSE ; CALL mpp_init2 ! eliminate land processors 286 306 ENDIF 307 #endif 287 308 ! 288 309 ! ! General initialization 310 ! CALL timing_init! Timing module 289 311 CALL phy_cst ! Physical constants 290 312 CALL eos_init ! Equation of state … … 482 504 USE trc_oce , ONLY: trc_oce_alloc 483 505 USE wrk_nemo , ONLY: wrk_alloc 506 USE exchmod , ONLY: exchmod_alloc 484 507 ! 485 508 INTEGER :: ierr … … 498 521 ierr = ierr + wrk_alloc(numout, lwp) ! workspace 499 522 ! 523 ierr = ierr + exchmod_alloc() ! New mpp msg framework 524 ! 500 525 IF( lk_mpp ) CALL mpp_sum( ierr ) 501 526 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'nemo_alloc : unable to allocate standard ocean arrays' ) … … 505 530 506 531 SUBROUTINE nemo_partition( num_pes ) 532 USE mapcomm_mod, ONLY: trimmed 507 533 !!---------------------------------------------------------------------- 508 534 !! *** ROUTINE nemo_partition *** … … 545 571 jpnij = jpni*jpnj 546 572 ! 573 574 ! Array that stores whether domain boundaries have been trimmed. Not used in 575 ! this case (regular domain decomp.) so set all to false. 576 ALLOCATE(trimmed(4,jpnij)) 577 trimmed(:,:) = .FALSE. 578 547 579 END SUBROUTINE nemo_partition 580 581 582 SUBROUTINE nemo_recursive_partition( num_pes ) 583 USE dom_oce, ONLY: ln_zco, ntopo 584 USE iom, ONLY: jpiglo, jpjglo, wp, jpdom_unknown, & 585 iom_open, iom_get, iom_close 586 USE mapcomm_mod, ONLY: ielb, ieub, pielb, pjelb, pieub, pjeub, & 587 iesub, jesub, jeub, ilbext, iubext, jubext, & 588 jlbext, pnactive, piesub, pjesub, jelb, pilbext, & 589 piubext, pjlbext, pjubext, LAND 590 USE partition_mod, ONLY: partition_rk, partition_mca_rk, imask, smooth_bathy 591 USE par_oce, ONLY: do_exchanges 592 #if defined key_mpp_mpi 593 USE mpi 594 #endif 595 !!---------------------------------------------------------------------- 596 !! *** ROUTINE nemo_recursive_partition *** 597 !! 598 !! ** Purpose : Work out a sensible factorisation of the number of 599 !! processors for the x and y dimensions. 600 !! ** Method : 601 !!---------------------------------------------------------------------- 602 IMPLICIT none 603 INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have 604 ! Local vars 605 INTEGER :: ierr ! Error flag 606 INTEGER :: inum ! temporary logical unit 607 INTEGER :: ii,jj,iproc ! Loop index 608 INTEGER :: jparray(2) ! Small array for gathering 609 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta ! temporary data workspace 610 !!---------------------------------------------------------------------- 611 612 ! Allocate masking array (stored in partition_mod) and workspace array 613 ! for this routine 614 ALLOCATE(imask(jpiglo,jpjglo), zdta(jpiglo,jpjglo), Stat=ierr) 615 IF(ierr /= 0)THEN 616 CALL ctl_stop('nemo_recursive_partition: failed to allocate workspace arrays') 617 RETURN 618 END IF 619 620 ! Factorise the number of MPI PEs to get jpi and jpj as usual 621 CALL nemo_partition(num_pes) 622 623 ! Generate a global mask... 624 !!$#if defined ARPDEBUG 625 !!$ IF(lwp)THEN 626 !!$ WRITE(*,*) 'ARPDBG: nemo_recursive_partition: generating mask...' 627 !!$ WRITE(*,*) 'ARPDBG: nemo_recursive_partition: jp{i,j}glo = ',jpiglo,jpjglo 628 !!$ END IF 629 !!$#endif 630 631 ! ARPDBG - this is the correct variable to check but the dom_nam section 632 ! of the namelist file hasn't been read in at this stage. 633 ! IF( ntopo == 1 )THEN 634 ! open the file 635 ierr = 0 636 !!$ IF ( ln_zco ) THEN 637 !!$ ! Setting ldstop prevents ctl_stop() from being called if the file 638 !!$ ! doesn't exist 639 !!$ CALL iom_open ( 'bathy_level.nc', inum, ldstop=.FALSE. ) ! Level bathymetry 640 !!$ IF(inum > 0)CALL iom_get ( inum, jpdom_unknown, 'Bathy_level', zdta, & 641 !!$ kstart=(/jpizoom,jpjzoom/), & 642 !!$ kcount=(/jpiglo,jpjglo/) ) 643 !!$ ELSE 644 CALL iom_open ( 'bathy_meter.nc', inum, ldstop=.FALSE. ) ! Meter bathy in case of partial steps 645 IF(inum > 0)CALL iom_get ( inum, jpdom_unknown, 'Bathymetry' , zdta, & 646 kstart=(/jpizoom,jpjzoom/), & 647 kcount=(/jpiglo,jpjglo/) ) 648 !!$ ENDIF 649 IF(inum > 0)THEN 650 CALL iom_close (inum) 651 ELSE 652 ! Flag that an error occurred when reading the file 653 ierr = 1 654 ENDIF 655 ! ELSE 656 ! ! Topography not read from file in this case 657 ! ierr = 1 658 ! END IF 659 660 ! If ln_sco defined then the bathymetry gets smoothed before the 661 ! simulation begins and that process can alter the coastlines 662 ! therefore we do it here too before calculating our mask. 663 ! IF(ln_sco) 664 CALL smooth_bathy(zdta) 665 666 ! land/sea mask (zero on land, 1 otherwise) over the global/zoom domain 667 imask(:,:)=1 668 IF(ierr == 1)THEN 669 ! Failed to read bathymetry so assume all ocean 670 WRITE(*,*) 'ARPDBG: nemo_recursive_partition: no bathymetry file so setting mask to unity' 671 672 ! Mess with otherwise uniform mask to get an irregular decomposition 673 ! for testing ARPDBG 674 CALL generate_fake_land(imask) 675 ELSE 676 ! Comment-out line below to achieve a regular partition 677 WHERE ( zdta(:,:) <= 1.0E-20 ) imask = LAND 678 END IF 679 680 ! Allocate partitioning arrays. 681 682 IF ( .not.allocated(pielb) ) THEN 683 ALLOCATE (pielb(num_pes), pieub(num_pes), piesub(num_pes), & 684 pilbext(num_pes), piubext(num_pes), & 685 pjelb(num_pes), pjeub(num_pes), pjesub(num_pes), & 686 pjlbext(num_pes), pjubext(num_pes), pnactive(num_pes), & 687 Stat = ierr) 688 IF(ierr /= 0)THEN 689 CALL ctl_stop('STOP', & 690 'nemo_recursive_partition: failed to allocate partitioning arrays') 691 RETURN 692 END IF 693 ENDIF 694 695 ! Now we can do recursive k-section partitioning 696 ! ARPDBG - BUG if limits on array below are set to anything other than 697 ! 1 and jp{i,j}glo then check for external boundaries in a few lines 698 ! time WILL FAIL! 699 ! CALL partition_rk ( imask, 1, jpiglo, 1, jpjglo, ierr ) 700 701 ! Multi-core aware version of recursive k-section partitioning 702 CALL partition_mca_rk ( imask, 1, jpiglo, 1, jpjglo, ierr ) 703 704 ! Check the error code from partitioning. 705 IF ( ierr /= 0 ) THEN 706 CALL ctl_stop('STOP','nemo_recursive_partition: Partitioning failed') 707 RETURN 708 ENDIF 709 710 ! Set the mask correctly now we've partitioned 711 !WHERE ( zdta(:,:) <= 0. ) imask = 0 712 713 ! ARPDBG Quick and dirty dump to stdout in gnuplot form 714 !!$ IF(narea == 1)THEN 715 !!$ OPEN(UNIT=998, FILE="imask.dat", & 716 !!$ STATUS='REPLACE', ACTION='WRITE', IOSTAT=jj) 717 !!$ IF( jj == 0 )THEN 718 !!$ WRITE (998,*) '# Depth map' 719 !!$ DO jj = 1, jpjglo, 1 720 !!$ DO ii = 1, jpiglo, 1 721 !!$ WRITE (998,*) ii, jj, zdta(ii,jj) ! imask(ii,jj) 722 !!$ END DO 723 !!$ WRITE (998,*) 724 !!$ END DO 725 !!$ CLOSE(998) 726 !!$ END IF 727 !!$ END IF 728 729 jpkm1 = jpk - 1 730 731 ! This chunk taken directly from original mpp_ini - not sure why nbondi 732 ! is reset? However, if it isn't reset then bad things happen in dommsk 733 ! so I'm doing what the original code does... 734 nperio = 0 735 nbondi = 0 736 IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN 737 IF( jpni == 1 )THEN 738 nbondi = 2 739 nperio = 1 740 END IF 741 END IF 742 743 #if defined ARPDEBUG 744 WRITE (*,FMT="(I4,' : ARPDBG: ielb, ieub, iesub = ',3I5)") narea-1,& 745 ielb, ieub, iesub 746 WRITE (*,FMT="(I4,' : ARPDBG: jelb, jeub, jesub = ',3I5)") narea-1,& 747 jelb, jeub, jesub 748 WRITE (*,FMT="(I4,' : ARPDBG: nldi, nlei, nlci = ',3I5)") narea-1, & 749 nldi, nlei, nlci 750 WRITE (*,FMT="(I4,' : ARPDBG: nldj, nlej, nlcj = ',3I5)") narea-1, & 751 nldj, nlej, nlcj 752 WRITE (*,FMT="(I4,' : ARPDBG: jpi, jpj = ',2I5)") narea-1, jpi, jpj 753 WRITE (*,FMT="(I4,' : ARPDBG: nimpp, njmpp = ',2I5)") narea-1, & 754 nimpp, njmpp 755 #endif 756 757 ! Debugging option - can turn off all halo exchanges by setting this to 758 ! false. 759 do_exchanges = .TRUE. 760 761 END SUBROUTINE nemo_recursive_partition 762 548 763 549 764 SUBROUTINE sqfact ( kn, kna, knb ) … … 565 780 566 781 fact_loop: DO kna=SQRT(REAL(kn)),1,-1 567 IF ( kn/kna*kna == kn ) THEN568 EXIT fact_loop569 ENDIF782 IF ( kn/kna*kna == kn ) THEN 783 EXIT fact_loop 784 ENDIF 570 785 END DO fact_loop 571 786 … … 577 792 END SUBROUTINE sqfact 578 793 794 795 SUBROUTINE generate_fake_land(imask) 796 !!---------------------------------------------------------------------- 797 !! Generate a fake land mass to test the decomposition code 798 !!---------------------------------------------------------------------- 799 USE par_oce, ONLY: jpiglo, jpjglo 800 USE partition_mod, ONLY: write_partition_map 801 IMPLICIT none 802 INTEGER, DIMENSION(jpiglo,jpjglo), INTENT(inout) :: imask 803 ! Locals 804 INTEGER :: ii, jj 805 INTEGER :: icentre, jcentre 806 INTEGER :: iwidth, iheight 807 INTEGER :: istart, istop 808 809 ! imask is zero on land points , unity on ocean points 810 iwidth = jpiglo/8 811 iheight = jpjglo/8 812 813 icentre = jpiglo/2 814 jcentre = jpjglo/2 815 816 istart = icentre - iwidth 817 istop = icentre + iwidth 818 DO jj = jcentre, jcentre - iheight, -1 819 imask(istart:istop,jj) = 0 820 istart = istart + 1 821 istop = istop - 1 822 END DO 823 istart = icentre - iwidth 824 istop = icentre + iwidth 825 DO jj = jcentre+1, jcentre + iheight, 1 826 imask(istart:istop,jj) = 0 827 istart = istart + 1 828 istop = istop - 1 829 END DO 830 831 ! Quick and dirty dump to stdout in gnuplot form 832 !!$ WRITE (*,*) 'GNUPLOT MAP' 833 !!$ DO jj = 1, jpjglo, 1 834 !!$ DO ii = 1, jpiglo, 1 835 !!$ WRITE (*,*) ii, jj, imask(ii,jj) 836 !!$ END DO 837 !!$ WRITE (*,*) 838 !!$ END DO 839 !!$ WRITE (*,*) 'END GNUPLOT MAP' 840 841 END SUBROUTINE generate_fake_land 842 579 843 !!====================================================================== 580 844 END MODULE nemogcm -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/par_oce.F90
r2715 r3432 81 81 !!--------------------------------------------------------------------- 82 82 # include "par_POMME_R025.h90" 83 #elif defined key_amm_12km 84 !!--------------------------------------------------------------------- 85 !! 'key_amm_12km': Atlantic Margin Model : AMM12km 86 !!--------------------------------------------------------------------- 87 # include "par_AMM_12km.h90" 83 88 #else 84 89 !!--------------------------------------------------------------------- … … 192 197 #endif 193 198 199 LOGICAL, PUBLIC :: do_exchanges ! ARPDBG for debugging only 200 ! If false then exch{r,s} 201 ! routines do nothing. 202 194 203 !!---------------------------------------------------------------------- 195 204 !! NEMO/OPA 3.3 , NEMO Consortium (2010) -
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/step.F90
r3187 r3432 37 37 #endif 38 38 USE asminc ! assimilation increments (tra_asm_inc, dyn_asm_inc routines) 39 39 USE timing, ONLY: timing_start, timing_stop, timing_reset, timing_disable 40 40 IMPLICIT NONE 41 41 PRIVATE … … 88 88 #endif 89 89 !! --------------------------------------------------------------------- 90 CALL timing_start('Step') 90 91 91 92 #if defined key_agrif … … 104 105 ! Update data, open boundaries, surface boundary condition (including sea-ice) 105 106 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 107 CALL timing_start('Boundaries') 106 108 IF( lk_dtatem ) CALL dta_tem( kstp ) ! update 3D temperature data 107 109 IF( lk_dtasal ) CALL dta_sal( kstp ) ! update 3D salinity data … … 110 112 IF( lk_obc ) CALL obc_rad( kstp ) ! compute phase velocities at open boundaries 111 113 IF( lk_bdy ) CALL bdy_dta_frs( kstp ) ! update dynamic and tracer data for FRS conditions (BDY) 114 ! ARP - no 'section' arg here so the time taken 115 ! in this region (which can involve IO) is 116 ! subtracted from the time 117 ! taken to do the whole step. 118 CALL timing_stop('Boundaries') 112 119 113 120 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 114 121 ! Ocean dynamics : ssh, wn, hdiv, rot ! 115 122 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 123 CALL timing_start('Ocean Dyn.') 116 124 CALL ssh_wzv( kstp ) ! after ssh & vertical velocity 125 CALL timing_stop('Ocean Dyn.','section') 117 126 118 127 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 119 128 ! Ocean physics update (ua, va, ta, sa used as workspace) 120 129 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 130 CALL timing_start('Ocean Phys.') 121 131 CALL bn2( tsb, rn2b ) ! before Brunt-Vaisala frequency 122 132 CALL bn2( tsn, rn2 ) ! now Brunt-Vaisala frequency 123 ! 124 ! VERTICAL PHYSICS 133 CALL timing_stop('Ocean Phys.','section') 134 ! 135 ! VERTICAL PHYSICS 136 CALL timing_start('Vert. Phys.') 125 137 CALL zdf_bfr( kstp ) ! bottom friction 126 138 … … 162 174 IF( lrst_oce .AND. lk_zdftke ) CALL tke_rst( kstp, 'WRITE' ) 163 175 IF( lrst_oce .AND. lk_zdfgls ) CALL gls_rst( kstp, 'WRITE' ) 176 177 CALL timing_stop('Vert. Phys.','section') 164 178 ! 165 179 ! LATERAL PHYSICS 180 ! 181 CALL timing_start('Lateral Phys.') 166 182 ! 167 183 IF( lk_ldfslp ) THEN ! slope of lateral mixing … … 178 194 IF( lk_traldf_eiv ) CALL ldf_eiv( kstp ) ! eddy induced velocity coefficient 179 195 #endif 196 CALL timing_stop('Lateral Phys.','section') 180 197 181 198 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 182 199 ! diagnostics and outputs (ua, va, ta, sa used as workspace) 183 200 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 201 CALL timing_start('Diagnostics') 184 202 IF( lk_floats ) CALL flo_stp( kstp ) ! drifting Floats 185 203 IF( lk_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20 degres isotherm depth) … … 188 206 IF( lk_diaar5 ) CALL dia_ar5( kstp ) ! ar5 diag 189 207 CALL dia_wri( kstp ) ! ocean model: outputs 208 CALL timing_stop('Diagnostics','section') 190 209 191 210 #if defined key_top … … 199 218 ! Active tracers (ua, va used as workspace) 200 219 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 220 CALL timing_start('Active tracers') 221 201 222 tsa(:,:,:,:) = 0.e0 ! set tracer trends to zero 202 223 … … 233 254 ENDIF 234 255 CALL tra_unswap ! udate T & S 3D arrays (to be suppressed) 256 CALL timing_stop('Active tracers','section') 235 257 236 258 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 237 259 ! Dynamics (ta, sa used as workspace) 238 260 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 261 CALL timing_start('Dynamics') 262 239 263 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 240 264 va(:,:,:) = 0.e0 … … 242 266 IF( ln_asmiau .AND. & 243 267 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 268 !CALL timing_start('dyn_adv') 244 269 CALL dyn_adv( kstp ) ! advection (vector or flux form) 270 !CALL timing_stop('dyn_adv','section') 271 272 !CALL timing_start('dyn_vor') 245 273 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 274 !CALL timing_stop('dyn_vor','section') 275 276 !CALL timing_start('dyn_ldf') 246 277 CALL dyn_ldf( kstp ) ! lateral mixing 278 !CALL timing_stop('dyn_ldf','section') 279 247 280 #if defined key_agrif 248 281 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge 249 282 #endif 283 !CALL timing_start('dyn_hpg') 250 284 CALL dyn_hpg( kstp ) &nb