- Timestamp:
- 2015-07-10T13:28:53+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r4624 r5581 8 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: 9 9 !! vvl option includes z_star and z_tilde coordinates 10 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 10 11 !!---------------------------------------------------------------------- 11 12 !! 'key_vvl' variable volume … … 44 45 45 46 !!* Namelist nam_vvl 46 LOGICAL , PUBLIC :: ln_vvl_zstar ! zstar vertical coordinate47 LOGICAL , PUBLIC :: ln_vvl_ztilde ! ztilde vertical coordinate48 LOGICAL , PUBLIC :: ln_vvl_layer ! level vertical coordinate49 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar ! ztilde vertical coordinate50 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor ! ztilde vertical coordinate51 LOGICAL , PUBLIC :: ln_vvl_kepe ! kinetic/potential energy transfer52 ! ! conservation: not used yet47 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate 48 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 49 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 50 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 51 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate 52 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 53 ! ! conservation: not used yet 53 54 REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient 54 55 REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days] 55 56 REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days] 56 57 REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation 57 LOGICAL , PUBLIC :: ln_vvl_dbg 58 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 58 59 59 60 !! * Module variables … … 125 126 INTEGER :: ji,jj,jk 126 127 INTEGER :: ii0, ii1, ij0, ij1 128 REAL(wp):: zcoef 127 129 !!---------------------------------------------------------------------- 128 130 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') … … 164 166 ! t- and w- points depth 165 167 ! ---------------------- 168 ! set the isf depth as it is in the initial step 166 169 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 167 170 fsdepw_n(:,:,1) = 0.0_wp … … 169 172 fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 170 173 fsdepw_b(:,:,1) = 0.0_wp 174 171 175 DO jk = 2, jpk 172 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 173 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 174 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 175 fsdept_b(:,:,jk) = fsdept_b(:,:,jk-1) + fse3w_b(:,:,jk) 176 fsdepw_b(:,:,jk) = fsdepw_b(:,:,jk-1) + fse3t_b(:,:,jk-1) 176 DO jj = 1,jpj 177 DO ji = 1,jpi 178 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 179 ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 180 ! 0.5 where jk = mikt 181 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 182 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 183 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 184 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 185 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 186 fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 187 fsdept_b(ji,jj,jk) = zcoef * ( fsdepw_b(ji,jj,jk ) + 0.5 * fse3w_b(ji,jj,jk)) & 188 & + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk)) 189 END DO 190 END DO 177 191 END DO 178 192 … … 185 199 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 186 200 END DO 187 hur_b(:,:) = umask (:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )188 hvr_b(:,:) = vmask (:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )201 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 202 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 189 203 190 204 ! Restoring frequencies for z_tilde coordinate … … 293 307 ! ! --------------------------------------------- ! 294 308 295 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )309 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 296 310 DO jk = 1, jpkm1 297 311 ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0) … … 313 327 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 314 328 END DO 315 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask (:,:,1) )329 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 316 330 317 331 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) … … 338 352 ELSE ! layer case 339 353 DO jk = 1, jpkm1 340 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) 354 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 341 355 END DO 342 356 END IF … … 386 400 ! d - thickness diffusion transport: boundary conditions 387 401 ! (stored for tracer advction and continuity equation) 388 CALL lbc_lnk( un_td , 'U' , -1. )389 CALL lbc_lnk( vn_td , 'V' , -1. )402 CALL lbc_lnk( un_td , 'U' , -1._wp) 403 CALL lbc_lnk( vn_td , 'V' , -1._wp) 390 404 391 405 ! 4 - Time stepping of baroclinic scale factors … … 398 412 z2dt = 2.0_wp * rdt 399 413 ENDIF 400 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )414 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 401 415 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 402 416 … … 453 467 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 454 468 END DO 455 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )469 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 456 470 DO jk = 1, jpkm1 457 471 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) … … 463 477 ! ! ---baroclinic part--------- ! 464 478 DO jk = 1, jpkm1 465 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) 479 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 466 480 END DO 467 481 ENDIF … … 531 545 END DO 532 546 ! ! Inverse of the local depth 533 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask (:,:,1) ) * umask(:,:,1)534 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask (:,:,1) ) * vmask(:,:,1)547 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 548 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 535 549 536 550 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) … … 569 583 INTEGER, INTENT( in ) :: kt ! time step 570 584 !! * Local declarations 571 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def572 INTEGER :: jk ! dummy loop indices585 INTEGER :: ji,jj,jk ! dummy loop indices 586 REAL(wp) :: zcoef 573 587 !!---------------------------------------------------------------------- 574 588 575 589 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp') 576 !577 CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def )578 590 ! 579 591 IF( kt == nit000 ) THEN … … 619 631 ! t- and w- points depth 620 632 ! ---------------------- 633 ! set the isf depth as it is in the initial step 621 634 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 622 635 fsdepw_n(:,:,1) = 0.0_wp 623 636 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 637 624 638 DO jk = 2, jpk 625 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 626 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 627 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 639 DO jj = 1,jpj 640 DO ji = 1,jpi 641 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 642 ! 1 for jk = mikt 643 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 644 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 645 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 646 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 647 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 648 END DO 649 END DO 628 650 END DO 651 629 652 ! Local depth and Inverse of the local depth of the water column at u- and v- points 630 653 ! ---------------------------------------------------------------------------------- … … 645 668 ! Write outputs 646 669 ! ============= 647 z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 648 CALL iom_put( "cellthc" , fse3t_n (:,:,:) ) 670 CALL iom_put( "e3t" , fse3t_n (:,:,:) ) 671 CALL iom_put( "e3u" , fse3u_n (:,:,:) ) 672 CALL iom_put( "e3v" , fse3v_n (:,:,:) ) 673 CALL iom_put( "e3w" , fse3w_n (:,:,:) ) 649 674 CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) ) 650 CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) ) 675 IF( iom_use("e3tdef") ) & 676 CALL iom_put( "e3tdef" , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 651 677 652 678 ! write restart file 653 679 ! ================== 654 680 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 655 !656 CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def )657 681 ! 658 682 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') … … 702 726 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 703 727 ! boundary conditions 704 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. )728 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 705 729 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 706 730 ! ! ------------------------------------- ! … … 720 744 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 721 745 ! boundary conditions 722 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. )746 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 723 747 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 724 748 ! ! ------------------------------------- ! … … 738 762 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 739 763 ! boundary conditions 740 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. )764 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 741 765 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 742 766 ! ! ------------------------------------- ! … … 808 832 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 809 833 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 810 id5 = iom_varid( numror, 'hdi f_lf', ldstop = .FALSE. )834 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 811 835 ! ! --------- ! 812 836 ! ! all cases ! … … 815 839 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 816 840 CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 841 ! needed to restart if land processor not computed 842 IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files' 843 WHERE ( tmask(:,:,:) == 0.0_wp ) 844 fse3t_n(:,:,:) = e3t_0(:,:,:) 845 fse3t_b(:,:,:) = e3t_0(:,:,:) 846 END WHERE 817 847 IF( neuler == 0 ) THEN 818 848 fse3t_b(:,:,:) = fse3t_n(:,:,:) 819 849 ENDIF 820 850 ELSE IF( id1 > 0 ) THEN 851 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files' 852 IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.' 853 IF(lwp) write(numout,*) 'neuler is forced to 0' 854 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 855 fse3t_n(:,:,:) = fse3t_b(:,:,:) 856 neuler = 0 857 ELSE IF( id2 > 0 ) THEN 821 858 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files' 822 859 IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.' 823 860 IF(lwp) write(numout,*) 'neuler is forced to 0' 861 CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 824 862 fse3t_b(:,:,:) = fse3t_n(:,:,:) 825 863 neuler = 0 … … 830 868 DO jk=1,jpk 831 869 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 832 & / ( ht_0(:,:) + 1._wp - tmask(:,:,1) ) * tmask(:,:,jk) &870 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 833 871 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 834 872 END DO … … 964 1002 965 1003 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 1004 IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 966 1005 967 1006 IF(lwp) THEN ! Print the choice … … 1000 1039 INTEGER :: ji, jj, jk ! dummy loop indices 1001 1040 INTEGER :: ij0, ij1, ii0, ii1 ! dummy loop indices 1041 INTEGER :: isrow ! index for ORCA1 starting row 1002 1042 !! acc 1003 1043 !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for … … 1083 1123 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration 1084 1124 ! ! ===================== 1085 ! 1086 ii0 = 281 ; ii1 = 282 ! Gibraltar Strait (e2u was modified) 1087 ij0 = 200 ; ij1 = 200 1125 ! This dirty section will be suppressed by simplification process: 1126 ! all this will come back in input files 1127 ! Currently these hard-wired indices relate to configuration with 1128 ! extend grid (jpjglo=332) 1129 ! which had a grid-size of 362x292. 1130 isrow = 332 - jpjglo 1131 ! 1132 ii0 = 282 ; ii1 = 283 ! Gibraltar Strait (e2u was modified) 1133 ij0 = 241 - isrow ; ij1 = 241 - isrow 1088 1134 DO jk = 1, jpkm1 1089 1135 DO jj = mj0(ij0), mj1(ij1) … … 1105 1151 END DO 1106 1152 ! 1107 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified)1108 ij0 = 2 08 ; ij1 = 2081153 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified) 1154 ij0 = 248 - isrow ; ij1 = 248 - isrow 1109 1155 DO jk = 1, jpkm1 1110 1156 DO jj = mj0(ij0), mj1(ij1) … … 1126 1172 END DO 1127 1173 ! 1128 ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified)1129 ij0 = 1 24 ; ij1 = 1251174 ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified) 1175 ij0 = 164 - isrow ; ij1 = 165 - isrow 1130 1176 DO jk = 1, jpkm1 1131 1177 DO jj = mj0(ij0), mj1(ij1) … … 1142 1188 END DO 1143 1189 ! 1144 ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on]1145 ij0 = 1 24 ; ij1 = 1251190 ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 1191 ij0 = 164 - isrow ; ij1 = 165 - isrow 1146 1192 DO jk = 1, jpkm1 1147 1193 DO jj = mj0(ij0), mj1(ij1) … … 1158 1204 END DO 1159 1205 ! 1160 ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified)1161 ij0 = 1 24 ; ij1 = 1251206 ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified) 1207 ij0 = 164 - isrow ; ij1 = 165 - isrow 1162 1208 DO jk = 1, jpkm1 1163 1209 DO jj = mj0(ij0), mj1(ij1) … … 1174 1220 END DO 1175 1221 ! 1176 ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified)1177 ij0 = 1 24 ; ij1 = 1251222 ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified) 1223 ij0 = 164 - isrow ; ij1 = 165 - isrow 1178 1224 DO jk = 1, jpkm1 1179 1225 DO jj = mj0(ij0), mj1(ij1) … … 1190 1236 END DO 1191 1237 ! 1192 ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified)1193 ij0 = 1 41 ; ij1 = 1421238 ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified) 1239 ij0 = 181 - isrow ; ij1 = 182 - isrow 1194 1240 DO jk = 1, jpkm1 1195 1241 DO jj = mj0(ij0), mj1(ij1) … … 1206 1252 END DO 1207 1253 ! 1208 ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified)1209 ij0 = 1 41 ; ij1 = 1421254 ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified) 1255 ij0 = 181 - isrow ; ij1 = 182 - isrow 1210 1256 DO jk = 1, jpkm1 1211 1257 DO jj = mj0(ij0), mj1(ij1)
Note: See TracChangeset
for help on using the changeset viewer.
