Changeset 13355
- Timestamp:
- 2020-07-30T12:12:41+02:00 (4 years ago)
- Location:
- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90
r13316 r13355 35 35 USE oce, ONLY: & ! active tracer variables 36 36 & tsn 37 USE zdfmxl, ONLY : & ! mixed layer depth 37 USE zdfmxl, ONLY : & ! mixed layer depth 38 38 #if defined key_karaml 39 39 & hmld_kara, & 40 40 & ln_kara, & 41 #endif 42 & hmld, & 41 #endif 42 & hmld, & 43 43 & hmlp, & 44 & hmlpt 44 & hmlpt 45 45 USE asmpar, ONLY: & ! assimilation parameters 46 46 & c_asmbkg, & … … 97 97 98 98 IMPLICIT NONE 99 PRIVATE 99 PRIVATE 100 100 101 101 PUBLIC asm_bgc_check_options ! called by asm_inc_init in asminc.F90 … … 304 304 305 305 ! Allocate and read increments 306 306 307 307 IF ( ln_slchltotinc ) THEN 308 308 ALLOCATE( slchltot_bkginc(jpi,jpj) ) 309 309 CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc ) 310 310 ENDIF 311 311 312 312 IF ( ln_slchldiainc ) THEN 313 313 ALLOCATE( slchldia_bkginc(jpi,jpj) ) 314 314 CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc ) 315 315 ENDIF 316 316 317 317 IF ( ln_slchlnoninc ) THEN 318 318 ALLOCATE( slchlnon_bkginc(jpi,jpj) ) 319 319 CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc ) 320 320 ENDIF 321 321 322 322 IF ( ln_schltotinc ) THEN 323 323 ALLOCATE( schltot_bkginc(jpi,jpj) ) 324 324 CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc ) 325 325 ENDIF 326 326 327 327 IF ( ln_slphytotinc ) THEN 328 328 ALLOCATE( slphytot_bkginc(jpi,jpj) ) 329 329 CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc ) 330 330 ENDIF 331 331 332 332 IF ( ln_slphydiainc ) THEN 333 333 ALLOCATE( slphydia_bkginc(jpi,jpj) ) 334 334 CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc ) 335 335 ENDIF 336 336 337 337 IF ( ln_slphynoninc ) THEN 338 338 ALLOCATE( slphynon_bkginc(jpi,jpj) ) … … 349 349 CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc ) 350 350 ENDIF 351 351 352 352 IF ( ln_plchltotinc ) THEN 353 353 ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) ) 354 354 CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc ) 355 355 ENDIF 356 356 357 357 IF ( ln_pchltotinc ) THEN 358 358 ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) ) 359 359 CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc ) 360 360 ENDIF 361 361 362 362 IF ( ln_pno3inc ) THEN 363 363 ALLOCATE( pno3_bkginc(jpi,jpj,jpk) ) 364 364 CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc ) 365 365 ENDIF 366 366 367 367 IF ( ln_psi4inc ) THEN 368 368 ALLOCATE( psi4_bkginc(jpi,jpj,jpk) ) 369 369 CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc ) 370 370 ENDIF 371 371 372 372 IF ( ln_pdicinc ) THEN 373 373 ALLOCATE( pdic_bkginc(jpi,jpj,jpk) ) 374 374 CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc ) 375 375 ENDIF 376 376 377 377 IF ( ln_palkinc ) THEN 378 378 ALLOCATE( palk_bkginc(jpi,jpj,jpk) ) 379 379 CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc ) 380 380 ENDIF 381 381 382 382 IF ( ln_pphinc ) THEN 383 383 ALLOCATE( pph_bkginc(jpi,jpj,jpk) ) 384 384 CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc ) 385 385 ENDIF 386 386 387 387 IF ( ln_po2inc ) THEN 388 388 ALLOCATE( po2_bkginc(jpi,jpj,jpk) ) … … 391 391 392 392 ! Allocate balancing increments 393 393 394 394 IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & 395 395 & ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & … … 402 402 #endif 403 403 ENDIF 404 404 405 405 IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN 406 406 #if defined key_top … … 457 457 ! Initialise 458 458 p_incs(:,:) = 0.0 459 459 460 460 ! read from file 461 461 CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 ) 462 462 463 463 ! Apply the masks 464 464 p_incs(:,:) = p_incs(:,:) * tmask(:,:,1) 465 465 466 466 ! Set missing increments to 0.0 rather than 1e+20 467 467 ! to allow for differences in masks … … 495 495 ! Initialise 496 496 p_incs(:,:,:) = 0.0 497 497 498 498 ! read from file 499 499 CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 ) 500 500 501 501 ! Apply the masks 502 502 p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:) 503 503 504 504 ! Set missing increments to 0.0 rather than 1e+20 505 505 ! to allow for differences in masks … … 558 558 cchl_p_bkg(:,:,:) = 0.0 559 559 #endif 560 560 561 561 !-------------------------------------------------------------------- 562 562 ! Read background variables for phytoplankton assimilation … … 578 578 CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) ) 579 579 #endif 580 580 581 581 IF ( ln_phytobal ) THEN 582 582 … … 628 628 629 629 CALL iom_close( inum ) 630 630 631 631 DO jt = 1, jptra 632 632 tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 633 633 END DO 634 634 635 635 ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN 636 636 … … 641 641 642 642 CALL iom_open( c_asmbkg, inum ) 643 643 644 644 #if defined key_hadocc 645 645 CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) … … 652 652 653 653 CALL iom_close( inum ) 654 654 655 655 DO jt = 1, jptra 656 656 tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 657 657 END DO 658 658 mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1) 659 659 660 660 ENDIF 661 661 #else … … 681 681 !! 682 682 !! ** Action : 683 !! 683 !! 684 684 !! References : 685 685 !! … … 695 695 REAL(wp) :: zdate ! Date 696 696 !!------------------------------------------------------------------------ 697 697 698 698 ! Set things up 699 699 zdate = REAL( ndastp ) … … 706 706 & TRIM( c_asmbal ) // ' at timestep = ', kt 707 707 708 ! Define the output file 708 ! Define the output file 709 709 CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib) 710 710 … … 812 812 & TRIM( c_asmbal ) // ' at timestep = ', kt 813 813 ENDIF 814 814 815 815 END SUBROUTINE asm_bgc_bal_wri 816 816 … … 893 893 !! ** Action : return non-log increments 894 894 !! 895 !! References : 895 !! References : 896 896 !!------------------------------------------------------------------------ 897 897 !! … … 970 970 !!------------------------------------------------------------------------ 971 971 !! *** ROUTINE phyto2d_asm_inc *** 972 !! 972 !! 973 973 !! ** Purpose : Apply the chlorophyll assimilation increments. 974 974 !! … … 977 977 !! Direct initialization or Incremental Analysis Updating. 978 978 !! 979 !! ** Action : 979 !! ** Action : 980 980 !!------------------------------------------------------------------------ 981 981 INTEGER, INTENT(IN) :: kt ! Current time step … … 1008 1008 REAL(wp), DIMENSION(jpi,jpj,1) :: zphyt_avg_bkg ! Local phyt_avg_bkg 1009 1009 !!------------------------------------------------------------------------ 1010 1010 1011 1011 IF ( kt <= nit000 ) THEN 1012 1012 … … 1014 1014 ! Remember that two sets of non-log increments should not be 1015 1015 ! expected to be in the same ratio as their log equivalents 1016 1016 1017 1017 ! Total chlorophyll 1018 1018 IF ( ln_slchltotinc ) THEN … … 1174 1174 1175 1175 IF(lwp) THEN 1176 WRITE(numout,*) 1176 WRITE(numout,*) 1177 1177 WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', & 1178 1178 & kt,' with IAU weight = ', pwgtiau(it) … … 1205 1205 ENDIF 1206 1206 1207 ELSEIF ( ll_asmdin ) THEN 1207 ELSEIF ( ll_asmdin ) THEN 1208 1208 1209 1209 !-------------------------------------------------------------------- 1210 1210 ! Direct Initialization 1211 1211 !-------------------------------------------------------------------- 1212 1212 1213 1213 IF ( kt == nitdin_r ) THEN 1214 1214 … … 1234 1234 END WHERE 1235 1235 #endif 1236 1236 1237 1237 ! Do not deallocate arrays - needed by asm_bgc_bal_wri 1238 1238 ! which is called at end of model run … … 1250 1250 !!------------------------------------------------------------------------ 1251 1251 !! *** ROUTINE phyto3d_asm_inc *** 1252 !! 1252 !! 1253 1253 !! ** Purpose : Apply the profile chlorophyll assimilation increments. 1254 1254 !! … … 1256 1256 !! Direct initialization or Incremental Analysis Updating. 1257 1257 !! 1258 !! ** Action : 1258 !! ** Action : 1259 1259 !!------------------------------------------------------------------------ 1260 1260 INTEGER, INTENT(IN) :: kt ! Current time step … … 1340 1340 1341 1341 IF(lwp) THEN 1342 WRITE(numout,*) 1342 WRITE(numout,*) 1343 1343 WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', & 1344 1344 & kt,' with IAU weight = ', pwgtiau(it) … … 1371 1371 ENDIF 1372 1372 1373 ELSEIF ( ll_asmdin ) THEN 1373 ELSEIF ( ll_asmdin ) THEN 1374 1374 1375 1375 !-------------------------------------------------------------------- 1376 1376 ! Direct Initialization 1377 1377 !-------------------------------------------------------------------- 1378 1378 1379 1379 IF ( kt == nitdin_r ) THEN 1380 1380 … … 1400 1400 END WHERE 1401 1401 #endif 1402 1402 1403 1403 ! Do not deallocate arrays - needed by asm_bgc_bal_wri 1404 1404 ! which is called at end of model run … … 1417 1417 !!------------------------------------------------------------------------ 1418 1418 !! *** ROUTINE pco2_asm_inc *** 1419 !! 1419 !! 1420 1420 !! ** Purpose : Apply the pco2/fco2 assimilation increments. 1421 1421 !! … … 1424 1424 !! Direct initialization or Incremental Analysis Updating. 1425 1425 !! 1426 !! ** Action : 1426 !! ** Action : 1427 1427 !!------------------------------------------------------------------------ 1428 1428 INTEGER, INTENT(IN) :: kt ! Current time step … … 1586 1586 jkmax = jpk-1 1587 1587 DO jk = jpk-1, 1, -1 1588 #if defined key_vvl 1588 1589 IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. & 1589 1590 & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN … … 1591 1592 jkmax = jk 1592 1593 ENDIF 1594 #else 1595 IF ( ( zmld(ji,jj) > gdepw_0(ji,jj,jk) ) .AND. & 1596 & ( zmld(ji,jj) <= gdepw_0(ji,jj,jk+1) ) ) THEN 1597 zmld(ji,jj) = gdepw_0(ji,jj,jk+1) 1598 jkmax = jk 1599 ENDIF 1600 #endif 1593 1601 END DO 1594 1602 ! … … 1617 1625 1618 1626 IF(lwp) THEN 1619 WRITE(numout,*) 1627 WRITE(numout,*) 1620 1628 IF ( ln_spco2inc ) THEN 1621 1629 WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', & … … 1654 1662 ENDIF 1655 1663 1656 ELSEIF ( ll_asmdin ) THEN 1664 ELSEIF ( ll_asmdin ) THEN 1657 1665 1658 1666 !-------------------------------------------------------------------- 1659 1667 ! Direct Initialization 1660 1668 !-------------------------------------------------------------------- 1661 1669 1662 1670 IF ( kt == nitdin_r ) THEN 1663 1671 … … 1683 1691 END WHERE 1684 1692 #endif 1685 1693 1686 1694 ! Do not deallocate arrays - needed by asm_bgc_bal_wri 1687 1695 ! which is called at end of model run … … 1700 1708 !!------------------------------------------------------------------------ 1701 1709 !! *** ROUTINE ph_asm_inc *** 1702 !! 1710 !! 1703 1711 !! ** Purpose : Apply the pH assimilation increments. 1704 1712 !! … … 1707 1715 !! Direct initialization or Incremental Analysis Updating. 1708 1716 !! 1709 !! ** Action : 1717 !! ** Action : 1710 1718 !!------------------------------------------------------------------------ 1711 1719 INTEGER, INTENT(IN) :: kt ! Current time step … … 1717 1725 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments 1718 1726 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments 1719 1727 1720 1728 REAL(wp) :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation 1721 1729 REAL(wp) :: DIC_IN, ALK_IN ! DIC/alk in pH calculation … … 1778 1786 sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) 1779 1787 ENDIF 1780 1788 1781 1789 ! Account for pCO2 balancing if required 1782 1790 IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN … … 1784 1792 alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk) 1785 1793 ENDIF 1786 1794 1787 1795 ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk 1788 1796 ! This requires three calls to the MOCSY carbonate package … … 1849 1857 ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic 1850 1858 ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk 1851 1859 1852 1860 ENDIF 1853 1861 1854 1862 END DO 1855 1863 END DO … … 1857 1865 1858 1866 ENDIF 1859 1867 1860 1868 IF ( ll_asmiau ) THEN 1861 1869 … … 1871 1879 1872 1880 IF(lwp) THEN 1873 WRITE(numout,*) 1881 WRITE(numout,*) 1874 1882 WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', & 1875 1883 & kt,' with IAU weight = ', pwgtiau(it) … … 1893 1901 ENDIF 1894 1902 1895 ELSEIF ( ll_asmdin ) THEN 1903 ELSEIF ( ll_asmdin ) THEN 1896 1904 1897 1905 !-------------------------------------------------------------------- 1898 1906 ! Direct Initialization 1899 1907 !-------------------------------------------------------------------- 1900 1908 1901 1909 IF ( kt == nitdin_r ) THEN 1902 1910 … … 1913 1921 trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 1914 1922 END WHERE 1915 1923 1916 1924 ! Do not deallocate arrays - needed by asm_bgc_bal_wri 1917 1925 ! which is called at end of model run … … 1919 1927 ! 1920 1928 ENDIF 1921 #endif 1929 #endif 1922 1930 ! 1923 1931 END SUBROUTINE ph_asm_inc … … 1930 1938 !!---------------------------------------------------------------------- 1931 1939 !! *** ROUTINE dyn_asm_inc *** 1932 !! 1940 !! 1933 1941 !! ** Purpose : Apply generic 3D biogeochemistry assimilation increments. 1934 1942 !! 1935 1943 !! ** Method : Direct initialization or Incremental Analysis Updating. 1936 1944 !! 1937 !! ** Action : 1945 !! ** Action : 1938 1946 !!---------------------------------------------------------------------- 1939 1947 INTEGER, INTENT(IN) :: kt ! Current time step … … 2082 2090 2083 2091 IF(lwp) THEN 2084 WRITE(numout,*) 2092 WRITE(numout,*) 2085 2093 WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', & 2086 2094 & kt,' with IAU weight = ', pwgtiau(it) … … 2170 2178 #endif 2171 2179 ENDIF 2172 2180 2173 2181 IF ( kt == nitiaufin_r ) THEN 2174 2182 IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc ) … … 2181 2189 ENDIF 2182 2190 2183 ELSEIF ( ll_asmdin ) THEN 2191 ELSEIF ( ll_asmdin ) THEN 2184 2192 2185 2193 !-------------------------------------------------------------------- 2186 2194 ! Direct Initialization 2187 2195 !-------------------------------------------------------------------- 2188 2196 2189 2197 IF ( kt == nitdin_r ) THEN 2190 2198 … … 2278 2286 #endif 2279 2287 ENDIF 2280 2288 2281 2289 IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc ) 2282 2290 IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc ) -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r8400 r13355 1 MODULE sshwzv 1 MODULE sshwzv 2 2 !!============================================================================== 3 3 !! *** MODULE sshwzv *** … … 5 5 !!============================================================================== 6 6 !! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code 7 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA 7 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA 8 8 !! - ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module … … 17 17 !!---------------------------------------------------------------------- 18 18 USE oce ! ocean dynamics and tracers variables 19 USE dom_oce ! ocean space and time domain variables 19 USE dom_oce ! ocean space and time domain variables 20 20 USE sbc_oce ! surface boundary condition: ocean 21 21 USE domvvl ! Variable volume … … 28 28 USE lib_mpp ! MPP library 29 29 USE bdy_oce 30 USE bdy_par 30 USE bdy_par 31 31 USE bdydyn2d ! bdy_ssh routine 32 32 #if defined key_agrif 33 33 USE agrif_opa_interp 34 34 #endif 35 #if defined key_asminc 35 #if defined key_asminc 36 36 USE asminc ! Assimilation increment 37 37 #endif … … 59 59 !!---------------------------------------------------------------------- 60 60 !! *** ROUTINE ssh_nxt *** 61 !! 61 !! 62 62 !! ** Purpose : compute the after ssh (ssha) 63 63 !! … … 73 73 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv 74 74 INTEGER, INTENT(in) :: kt ! time step 75 ! 75 ! 76 76 INTEGER :: jk ! dummy loop indices 77 77 REAL(wp) :: z2dt, z1_rau0 ! local scalars … … 80 80 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 81 81 ! 82 CALL wrk_alloc( jpi, jpj, zhdiv ) 82 CALL wrk_alloc( jpi, jpj, zhdiv ) 83 83 ! 84 84 IF( kt == nit000 ) THEN … … 102 102 #if defined key_vvl 103 103 ! Don't directly adjust ssh but change hdivn at all levels instead 104 ! In trasbc also add in the heat and salt content associated with these changes at each level 105 DO jk = 1, jpkm1 106 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk) 104 ! In trasbc also add in the heat and salt content associated with these changes at each level 105 DO jk = 1, jpkm1 106 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk) 107 107 END DO 108 ENDIF 109 #endif 108 #endif 109 ENDIF 110 110 #endif 111 111 … … 121 121 ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to 122 122 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 123 ! 123 ! 124 124 z1_rau0 = 0.5_wp * r1_rau0 125 125 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) … … 146 146 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 147 147 ! 148 CALL wrk_dealloc( jpi, jpj, zhdiv ) 148 CALL wrk_dealloc( jpi, jpj, zhdiv ) 149 149 ! 150 150 IF( nn_timing == 1 ) CALL timing_stop('ssh_nxt') … … 152 152 END SUBROUTINE ssh_nxt 153 153 154 154 155 155 SUBROUTINE wzv( kt ) 156 156 !!---------------------------------------------------------------------- 157 157 !! *** ROUTINE wzv *** 158 !! 158 !! 159 159 !! ** Purpose : compute the now vertical velocity 160 160 !! 161 !! ** Method : - Using the incompressibility hypothesis, the vertical 162 !! velocity is computed by integrating the horizontal divergence 161 !! ** Method : - Using the incompressibility hypothesis, the vertical 162 !! velocity is computed by integrating the horizontal divergence 163 163 !! from the bottom to the surface minus the scale factor evolution. 164 164 !! The boundary conditions are w=0 at the bottom (no flux) and. … … 176 176 REAL(wp) :: z1_2dt ! local scalars 177 177 !!---------------------------------------------------------------------- 178 178 179 179 IF( nn_timing == 1 ) CALL timing_start('wzv') 180 180 ! … … 195 195 ! 196 196 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases 197 CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 197 CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 198 198 ! 199 199 DO jk = 1, jpkm1 … … 215 215 END DO 216 216 ! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 217 CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 217 CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 218 218 ELSE ! z_star and linear free surface cases 219 219 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence … … 241 241 !! *** ROUTINE ssh_nxt *** 242 242 !! 243 !! ** Purpose : achieve the sea surface height time stepping by 243 !! ** Purpose : achieve the sea surface height time stepping by 244 244 !! applying Asselin time filter and swapping the arrays 245 !! ssha already computed in ssh_nxt 245 !! ssha already computed in ssh_nxt 246 246 !! 247 247 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
r11442 r13355 21 21 USE lib_mpp ! MPP library 22 22 USE wrk_nemo ! work arrays 23 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 23 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 24 24 USE stopack 25 25 … … 31 31 32 32 INTEGER :: albd_init = 0 !: control flag for initialization 33 33 34 34 REAL(wp) :: rmue = 0.40 ! cosine of local solar altitude 35 35 REAL(wp) :: ralb_oce = 0.066 ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) … … 37 37 REAL(wp) :: c2 = 0.10 ! " " 38 38 REAL(wp) :: rcloud = 0.06 ! cloud effect on albedo (only-for nn_ice_alb=0) 39 39 40 40 ! !!* namelist namsbc_alb 41 41 INTEGER :: nn_ice_alb … … 52 52 !!---------------------------------------------------------------------- 53 53 !! *** ROUTINE albedo_ice *** 54 !! 55 !! ** Purpose : Computation of the albedo of the snow/ice system 56 !! 54 !! 55 !! ** Purpose : Computation of the albedo of the snow/ice system 56 !! 57 57 !! ** Method : Two schemes are available (from namelist parameter nn_ice_alb) 58 58 !! 0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies … … 73 73 !! ** Note : The parameterization from Shine & Henderson-Sellers presents several misconstructions: 74 74 !! 1) ice albedo when ice thick. tends to 0 is different than ocean albedo 75 !! 2) for small ice thick. covered with some snow (<3cm?), albedo is larger 75 !! 2) for small ice thick. covered with some snow (<3cm?), albedo is larger 76 76 !! under melting conditions than under freezing conditions 77 !! 3) the evolution of ice albedo as a function of ice thickness shows 77 !! 3) the evolution of ice albedo as a function of ice thickness shows 78 78 !! 3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic 79 79 !! 80 80 !! References : Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250. 81 81 !! Brandt et al. 2005, J. Climate, vol 18 82 !! Grenfell & Perovich 2004, JGR, vol 109 82 !! Grenfell & Perovich 2004, JGR, vol 109 83 83 !!---------------------------------------------------------------------- 84 84 REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: pt_ice ! ice surface temperature (Kelvin) … … 97 97 98 98 ijpl = SIZE( pt_ice, 3 ) ! number of ice categories 99 99 100 100 CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it ) 101 101 102 IF( albd_init == 0 ) CALL albedo_init ! initialization 103 104 102 IF( albd_init == 0 ) CALL albedo_init ! initialization 103 104 105 105 SELECT CASE ( nn_ice_alb ) 106 106 … … 109 109 !------------------------------------------ 110 110 CASE( 0 ) 111 111 112 112 ralb_sf = 0.80 ! dry snow 113 113 ralb_sm = 0.65 ! melting snow 114 114 ralb_if = 0.72 ! bare frozen ice 115 ralb_im = rn_albice ! bare puddled ice 116 115 ralb_im = rn_albice ! bare puddled ice 116 117 117 ! Computation of ice albedo (free of snow) 118 118 WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice ) ; zalb(:,:,:) = ralb_im 119 119 ELSE WHERE ; zalb(:,:,:) = ralb_if 120 120 END WHERE 121 121 122 122 WHERE ( 1.5 < ph_ice ) ; zalb_it = zalb 123 123 ELSE WHERE( 1.0 < ph_ice .AND. ph_ice <= 1.5 ) ; zalb_it = 0.472 + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) … … 127 127 ELSE WHERE ; zalb_it = 0.1 + 3.6 * ph_ice 128 128 END WHERE 129 129 130 130 DO jl = 1, ijpl 131 131 DO jj = 1, jpj … … 133 133 ! freezing snow 134 134 ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2 135 ! ! freezing snow 135 ! ! freezing snow 136 136 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 137 137 zalb_sf = ( 1._wp - zswitch ) * ( zalb_it(ji,jj,jl) & 138 138 & + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1 ) & 139 & + zswitch * ralb_sf 139 & + zswitch * ralb_sf 140 140 141 141 ! melting snow … … 143 143 zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) ) 144 144 zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 ) & 145 & + zswitch * ralb_sm 145 & + zswitch * ralb_sm 146 146 ! 147 147 ! snow albedo 148 zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) ) 148 zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) ) 149 149 zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 150 150 151 151 ! Ice/snow albedo 152 152 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) … … 155 155 END DO 156 156 END DO 157 157 158 #if defined key_traldf_c2d || key_traldf_c3d 158 159 IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 159 160 & CALL spp_gen( 1, pa_ice_cs(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl ) 160 161 #else 162 IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 163 & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// & 164 'key_traldf_c2d or key_traldf_c3d') 165 #endif 161 166 END DO 162 167 … … 166 171 ! New parameterization (2016) 167 172 !------------------------------------------ 168 CASE( 1 ) 173 CASE( 1 ) 169 174 170 175 ralb_im = rn_albice ! bare puddled ice … … 181 186 ! ralb_sm = 0.82 ! melting snow 182 187 ! ralb_if = 0.54 ! bare frozen ice 183 ! 188 ! 184 189 ! Computation of ice albedo (free of snow) 185 z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 190 z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 186 191 z1_c2 = 1. / 0.05 187 192 WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice ) ; zalb = ralb_im 188 193 ELSE WHERE ; zalb = ralb_if 189 194 END WHERE 190 195 191 196 WHERE ( 1.5 < ph_ice ) ; zalb_it = zalb 192 197 ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 ) ; zalb_it = zalb + ( 0.18 - zalb ) * z1_c1 * & … … 205 210 206 211 ! snow albedo 207 zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) ) 212 zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) ) 208 213 zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 209 214 210 ! Ice/snow albedo 215 ! Ice/snow albedo 211 216 zswitch = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 212 217 pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch * zalb_it(ji,jj,jl) … … 214 219 END DO 215 220 END DO 216 221 222 #if defined key_traldf_c2d || key_traldf_c3d 217 223 IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 218 224 & CALL spp_gen( 1, pa_ice_os(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl ) 219 225 #else 226 IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 227 & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// & 228 'key_traldf_c2d or key_traldf_c3d') 229 #endif 220 230 END DO 221 231 ! Effect of the clouds (2d order polynomial) 222 pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 232 pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 223 233 224 234 END SELECT 225 235 226 236 CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it ) 227 237 ! … … 232 242 !!---------------------------------------------------------------------- 233 243 !! *** ROUTINE albedo_oce *** 234 !! 244 !! 235 245 !! ** Purpose : Computation of the albedo of the ocean 236 246 !!---------------------------------------------------------------------- … … 238 248 REAL(wp), DIMENSION(:,:), INTENT(out) :: pa_oce_cs ! albedo of ocean under clear sky 239 249 !! 240 REAL(wp) :: zcoef 250 REAL(wp) :: zcoef 241 251 !!---------------------------------------------------------------------- 242 252 ! 243 253 zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 ) ! Parameterization of Briegled and Ramanathan, 1982 244 pa_oce_cs(:,:) = zcoef 254 pa_oce_cs(:,:) = zcoef 245 255 pa_oce_os(:,:) = 0.06 ! Parameterization of Kondratyev, 1969 and Payne, 1972 246 256 ! … … 257 267 !!---------------------------------------------------------------------- 258 268 INTEGER :: ios ! Local integer output status for namelist read 259 NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 269 NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 260 270 !!---------------------------------------------------------------------- 261 271 ! -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r11442 r13355 161 161 ! ! ====================== ! 162 162 ! 163 rn_sfac = 1._wp ! Default to one if missing from namelist 163 rn_sfac = 1._wp ! Default to one if missing from namelist 164 164 REWIND( numnam_ref ) ! Namelist namsbc_core in reference namelist : CORE bulk parameters 165 165 READ ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901) … … 205 205 ENDIF 206 206 207 #if defined key_traldf_c2d || key_traldf_c3d 207 208 IF( ln_stopack .AND. nn_spp_relw > 0 ) THEN 208 209 rn_vfac0(:,:) = rn_vfac 209 210 CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw ) 210 211 ENDIF 212 #else 213 IF ( ln_stopack .AND. nn_spp_relw > 0 ) & 214 & CALL ctl_stop( 'sbc_blk_core: parameter perturbation will only work with '// & 215 'key_traldf_c2d or key_traldf_c3d') 216 #endif 211 217 212 218 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step … … 217 223 #if defined key_cice 218 224 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 219 qlw_ice(:,:,1) = sf(jp_qlw)%fnow(:,:,1) 225 qlw_ice(:,:,1) = sf(jp_qlw)%fnow(:,:,1) 220 226 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 221 227 ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 222 228 ENDIF 223 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 229 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 224 230 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 225 231 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac … … 231 237 ! 232 238 END SUBROUTINE sbc_blk_core 233 234 239 240 235 241 SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv ) 236 242 !!--------------------------------------------------------------------- … … 242 248 !! ** Method : CORE bulk formulea for the ocean using atmospheric 243 249 !! fields read in sbc_read 244 !! 250 !! 245 251 !! ** Outputs : - utau : i-component of the stress at U-point (N/m2) 246 252 !! - vtau : j-component of the stress at V-point (N/m2) … … 280 286 ! local scalars ( place there for vector optimisation purposes) 281 287 zcoef_qsatw = 0.98 * 640380. / rhoa 282 288 283 289 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 284 290 … … 288 294 289 295 ! ... components ( U10m - U_oce ) at T-point (unmasked) 290 zwnd_i(:,:) = 0.e0 296 zwnd_i(:,:) = 0.e0 291 297 zwnd_j(:,:) = 0.e0 292 298 #if defined key_cyclone … … 332 338 CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm, & 333 339 & Cd, Ch, Ce, zt_zu, zq_zu ) 334 340 335 341 ! ... tau module, i and j component 336 342 DO jj = 1, jpj … … 363 369 CALL lbc_lnk( vtau(:,:), 'V', -1. ) 364 370 365 371 366 372 ! Turbulent fluxes over ocean 367 373 ! ----------------------------- … … 388 394 CALL prt_ctl( tab2d_1=zst , clinfo1=' blk_oce_core: zst : ') 389 395 ENDIF 390 396 391 397 ! ----------------------------------------------------------------------------- ! 392 398 ! III Total FLUXES ! … … 396 402 & - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1) 397 403 ! 398 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar 404 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar 399 405 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus & ! remove latent melting heat for solid precip 400 406 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST … … 437 443 ! 438 444 END SUBROUTINE blk_oce_core 439 440 445 446 441 447 #if defined key_lim2 || defined key_lim3 442 448 SUBROUTINE blk_ice_core_tau … … 531 537 532 538 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core_tau') 533 539 534 540 END SUBROUTINE blk_ice_core_tau 535 541 … … 544 550 !! between atmosphere and sea-ice using CORE bulk 545 551 !! formulea, ice variables and read atmmospheric fields. 546 !! 552 !! 547 553 !! caution : the net upward water flux has with mm/day unit 548 554 !!--------------------------------------------------------------------- … … 564 570 IF( nn_timing == 1 ) CALL timing_start('blk_ice_core_flx') 565 571 ! 566 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 572 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 567 573 568 574 ! local scalars ( place there for vector optimisation purposes) … … 587 593 z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 588 594 ! lw sensitivity 589 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 595 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 590 596 591 597 ! ----------------------------! … … 597 603 z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 598 604 ! Latent Heat 599 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * Cice * wndm_ice(ji,jj) & 605 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * Cice * wndm_ice(ji,jj) & 600 606 & * ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1) ) ) 601 607 ! Latent heat sensitivity for ice (Dqla/Dt) … … 628 634 629 635 #if defined key_lim3 630 CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 636 CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 631 637 632 638 ! --- evaporation --- ! … … 638 644 ! --- evaporation minus precipitation --- ! 639 645 zsnw(:,:) = 0._wp 640 CALL lim_thd_snwblow( pfrld, zsnw ) ! snow distribution over ice after wind blowing 646 CALL lim_thd_snwblow( pfrld, zsnw ) ! snow distribution over ice after wind blowing 641 647 emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 642 648 emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw … … 661 667 DO jl = 1, jpl 662 668 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 663 ! But we do not have Tice => consider it at 0 degC => evap=0 669 ! But we do not have Tice => consider it at 0 degC => evap=0 664 670 END DO 665 671 666 CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 672 CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 667 673 #endif 668 674 … … 688 694 ! 689 695 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core_flx') 690 696 691 697 END SUBROUTINE blk_ice_core_flx 692 698 #endif … … 701 707 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 702 708 !! 703 !! ** Method : Monin Obukhov Similarity Theory 709 !! ** Method : Monin Obukhov Similarity Theory 704 710 !! + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 705 711 !! … … 711 717 !! - better first guess of stability by checking air-sea difference of virtual temperature 712 718 !! rather than temperature difference only... 713 !! - added function "cd_neutral_10m" that uses the improved parametrization of 719 !! - added function "cd_neutral_10m" that uses the improved parametrization of 714 720 !! Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 715 721 !! - using code-wide physical constants defined into "phycst.mod" rather than redifining them … … 746 752 747 753 IF( nn_timing == 1 ) CALL timing_start('turb_core_2z') 748 754 749 755 CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 750 756 CALL wrk_alloc( jpi,jpj, zeta_u, stab ) … … 758 764 U_zu = MAX( 0.5 , dU ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 759 765 760 !! First guess of stability: 766 !! First guess of stability: 761 767 ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt 762 768 stab = 0.5 + sign(0.5,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable … … 772 778 Ce_n10 = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 773 779 Ch_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 774 780 775 781 !! Initializing transf. coeff. with their first guess neutral equivalents : 776 782 Cd = ztmp0 ; Ce = Ce_n10 ; Ch = Ch_n10 ; sqrt_Cd = sqrt_Cd_n10 … … 783 789 ! 784 790 ztmp1 = T_zu - sst ! Updating air/sea differences 785 ztmp2 = q_zu - q_sat 791 ztmp2 = q_zu - q_sat 786 792 787 793 ! Updating turbulent scales : (L&Y 2004 eq. (7)) 788 794 ztmp1 = Ch/sqrt_Cd*ztmp1 ! theta* 789 795 ztmp2 = Ce/sqrt_Cd*ztmp2 ! q* 790 796 791 797 ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu 792 798 793 799 ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 794 ztmp0 = (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 800 ztmp0 = (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 795 801 ! ( Cd*U_zu*U_zu is U*^2 at zu) 796 802 … … 799 805 zpsi_h_u = psi_h( zeta_u ) 800 806 zpsi_m_u = psi_m( zeta_u ) 801 807 802 808 !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c)) 803 809 IF ( .NOT. l_zt_equal_zu ) THEN … … 808 814 q_zu = max(0., q_zu) 809 815 END IF 810 816 811 817 IF( ln_cdgw ) THEN ! surface wave case 812 sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 818 sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 813 819 Cd = sqrt_Cd * sqrt_Cd 814 820 ELSE … … 820 826 ztmp0 = cd_neutral_10m(ztmp0) ! Cd_n10 821 827 sqrt_Cd_n10 = sqrt(ztmp0) 822 828 823 829 Ce_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) 824 830 stab = 0.5 + sign(0.5,zeta_u) ! update stability … … 827 833 !! Update of transfer coefficients: 828 834 ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u) ! L&Y 2004 eq. (10a) 829 Cd = ztmp0 / ( ztmp1*ztmp1 ) 835 Cd = ztmp0 / ( ztmp1*ztmp1 ) 830 836 sqrt_Cd = SQRT( Cd ) 831 837 ENDIF … … 833 839 ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 834 840 ztmp2 = sqrt_Cd / sqrt_Cd_n10 835 ztmp1 = 1. + Ch_n10*ztmp0 841 ztmp1 = 1. + Ch_n10*ztmp0 836 842 Ch = Ch_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 837 843 ! 838 ztmp1 = 1. + Ce_n10*ztmp0 844 ztmp1 = 1. + Ce_n10*ztmp0 839 845 Ce = Ce_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 840 846 ! … … 854 860 FUNCTION cd_neutral_10m( zw10 ) 855 861 !!---------------------------------------------------------------------- 856 !! Estimate of the neutral drag coefficient at 10m as a function 862 !! Estimate of the neutral drag coefficient at 10m as a function 857 863 !! of neutral wind speed at 10m 858 864 !! … … 860 866 !! 861 867 !! Author: L. Brodeau, june 2014 862 !!---------------------------------------------------------------------- 868 !!---------------------------------------------------------------------- 863 869 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zw10 ! scalar wind speed at 10m (m/s) 864 870 REAL(wp), DIMENSION(jpi,jpj) :: cd_neutral_10m 865 871 ! 866 872 REAL(wp), DIMENSION(:,:), POINTER :: rgt33 867 !!---------------------------------------------------------------------- 873 !!---------------------------------------------------------------------- 868 874 ! 869 875 CALL wrk_alloc( jpi,jpj, rgt33 ) 870 876 ! 871 877 !! When wind speed > 33 m/s => Cyclone conditions => special treatment 872 rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) ) ! If zw10 < 33. => 0, else => 1 878 rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) ) ! If zw10 < 33. => 0, else => 1 873 879 cd_neutral_10m = 1.e-3 * ( & 874 880 & (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
r11442 r13355 24 24 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 25 25 USE timing ! Timing 26 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 26 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 27 27 USE stopack 28 28 USE wrk_nemo ! Memory Allocation … … 42 42 REAL(wp) :: rn_dqdt ! restoring factor on SST and SSS 43 43 REAL(wp) :: rn_deds ! restoring factor on SST and SSS 44 LOGICAL :: ln_sssr_bnd ! flag to bound erp term 44 LOGICAL :: ln_sssr_bnd ! flag to bound erp term 45 45 REAL(wp) :: rn_sssr_bnd ! ABS(Max./Min.) value of erp term [mm/day] 46 46 … … 101 101 rn_dqdt_s(:,:) = rn_dqdt 102 102 103 IF( ln_stopack .AND. nn_spp_dqdt > 0 ) & 104 & CALL spp_gen( kt, rn_dqdt_s, nn_spp_dqdt, rn_dqdt_sd, jk_spp_dqdt ) 103 #if defined key_traldf_c2d || key_traldf_c3d 104 IF( ln_stopack .AND. nn_spp_dqdt > 0 ) & 105 & CALL spp_gen( kt, rn_dqdt_s, nn_spp_dqdt, rn_dqdt_sd, jk_spp_dqdt ) 106 #else 107 IF ( ln_stopack .AND. nn_spp_dqdt > 0 ) & 108 & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// & 109 'key_traldf_c2d or key_traldf_c3d') 110 #endif 111 105 112 DO jj = 1, jpj 106 113 DO ji = 1, jpi … … 117 124 CALL wrk_alloc( jpi, jpj, zsrp) 118 125 zsrp(:,:) = rn_deds 126 #if defined key_traldf_c2d || key_traldf_c3d 119 127 IF( ln_stopack .AND. nn_spp_dedt > 0 ) & 120 128 & CALL spp_gen(kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds ) 129 #else 130 IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 131 & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// & 132 'key_traldf_c2d or key_traldf_c3d') 133 #endif 134 135 121 136 !CDIR COLLAPSE 122 137 DO jj = 1, jpj 123 138 DO ji = 1, jpi 124 139 zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) ) & ! No damping in vicinity of river mouths 125 & * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 140 & * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 126 141 sfx(ji,jj) = sfx(ji,jj) + zerp ! salt flux 127 142 erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only) … … 134 149 CALL wrk_alloc( jpi, jpj, zsrp) 135 150 zsrp(:,:) = rn_deds 151 #if defined key_traldf_c2d || key_traldf_c3d 136 152 IF( ln_stopack .AND. nn_spp_dedt > 0 ) & 137 153 & CALL spp_gen( kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds ) 138 zerp_bnd = rn_sssr_bnd / rday ! - - 154 #else 155 IF ( ln_stopack .AND. nn_spp_dedt > 0 ) & 156 & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// & 157 'key_traldf_c2d or key_traldf_c3d') 158 #endif 159 zerp_bnd = rn_sssr_bnd / rday ! - - 139 160 !CDIR COLLAPSE 140 161 DO jj = 1, jpj 141 DO ji = 1, jpi 162 DO ji = 1, jpi 142 163 zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) ) & ! No damping in vicinity of river mouths 143 164 & * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) & … … 161 182 END SUBROUTINE sbc_ssr 162 183 163 184 164 185 SUBROUTINE sbc_ssr_init 165 186 !!--------------------------------------------------------------------- … … 184 205 !!---------------------------------------------------------------------- 185 206 ! 186 187 REWIND( numnam_ref ) ! Namelist namsbc_ssr in reference namelist : 207 208 REWIND( numnam_ref ) ! Namelist namsbc_ssr in reference namelist : 188 209 READ ( numnam_ref, namsbc_ssr, IOSTAT = ios, ERR = 901) 189 210 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_ssr in reference namelist', lwp ) … … 240 261 ENDIF 241 262 ! 242 ! !* Initialize qrp and erp if no restoring 263 ! !* Initialize qrp and erp if no restoring 243 264 IF( nn_sstr /= 1 ) qrp(:,:) = 0._wp 244 265 IF( nn_sssr /= 1 .OR. nn_sssr /= 2 ) erp(:,:) = 0._wp 245 266 ! 246 267 END SUBROUTINE sbc_ssr_init 247 268 248 269 !!====================================================================== 249 270 END MODULE sbcssr -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/STO/stopack.F90
r12102 r13355 46 46 #define numnam_ref numnam 47 47 #define numnam_cfg numnam 48 #define lwm lwp 48 #define lwm lwp 49 49 #define numond numout 50 50 51 #define wmask tmask 51 #define wmask tmask 52 52 53 53 #endif … … 63 63 !! (SPP, SKEB and sea-ice) 64 64 !!---------------------------------------------------------------------- 65 !! 65 !! 66 66 !! stopack : Generate stochastic physics perturbations 67 !! 67 !! 68 68 !! Method 69 69 !! ====== … … 71 71 !! - SPPT (Stochastically perturbed parameterization 72 72 !! tendencies )scheme for user-selected trends for 73 !! tracers, momentum and sea-ice 73 !! tracers, momentum and sea-ice 74 74 !! - SPP (Schastically perturbed parameters) scheme 75 75 !! for some (namelist) parameters … … 77 77 !! backscatter energy dissipated numerically or 78 78 !! through deep convection. 79 !! 80 !! 79 !! 80 !! 81 81 !! Acknowledgements: C3S funded ERGO project 82 !! 82 !! 83 83 !!---------------------------------------------------------------------- 84 84 USE par_kind 85 85 USE timing ! Timing 86 86 USE oce ! ocean dynamics and tracers variables 87 USE dom_oce ! ocean space and time domain variables 87 USE dom_oce ! ocean space and time domain variables 88 88 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 89 89 USE in_out_manager ! I/O manager … … 104 104 USE wrk_nemo 105 105 USE diaptr 106 USE zdf_oce 106 USE zdf_oce 107 107 USE phycst 108 108 … … 113 113 PUBLIC tra_sppt_collect 114 114 PUBLIC dyn_sppt_collect 115 PUBLIC tra_sppt_apply 116 PUBLIC dyn_sppt_apply 115 PUBLIC tra_sppt_apply 116 PUBLIC dyn_sppt_apply 117 117 PUBLIC stopack_rst 118 118 PUBLIC stopack_init … … 140 140 REAL(wp), SAVE :: rn_spp_tau, rn_spp_stdev 141 141 INTEGER :: skeb_filter_pass, spp_filter_pass 142 142 143 143 ! SPPT Logical switches for individual tendencies 144 144 LOGICAL :: ln_sppt_taumap, ln_stopack_restart, ln_distcoast, & 145 145 ln_sppt_traxad, ln_sppt_trayad, ln_sppt_trazad, ln_sppt_trasad, ln_sppt_traldf, & 146 146 ln_sppt_trazdf, ln_sppt_trazdfp,ln_sppt_traevd, ln_sppt_trabbc, ln_sppt_trabbl, & 147 ln_sppt_tranpc, ln_sppt_tradmp, ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf 147 ln_sppt_tranpc, ln_sppt_tradmp, ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf 148 148 LOGICAL :: & 149 149 ln_sppt_dynhpg, ln_sppt_dynspg, ln_sppt_dynkeg, ln_sppt_dynrvo, ln_sppt_dynpvo, ln_sppt_dynzad,& … … 181 181 INTEGER, PARAMETER, PUBLIC :: jk_spp_qsi0 = 8 182 182 INTEGER, PARAMETER, PUBLIC :: jk_spp_bfr = 9 183 INTEGER, PARAMETER, PUBLIC :: jk_spp_aevd = 10 183 INTEGER, PARAMETER, PUBLIC :: jk_spp_aevd = 10 184 184 INTEGER, PARAMETER, PUBLIC :: jk_spp_avt = 11 185 185 INTEGER, PARAMETER, PUBLIC :: jk_spp_avm = 12 … … 219 219 INTEGER, SAVE :: numrep = 602 220 220 INTEGER, SAVE :: lkt 221 221 222 222 ! Randome generator seed 223 223 INTEGER, SAVE :: nn_stopack_seed(4) … … 275 275 !!---------------------------------------------------------------------- 276 276 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 277 !! $Id: bdytra.F90 4292 2013-11-20 16:28:04Z cetlod $ 277 !! $Id: bdytra.F90 4292 2013-11-20 16:28:04Z cetlod $ 278 278 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 279 279 !!---------------------------------------------------------------------- … … 289 289 !! 290 290 !! ** Purpose : Collect tracer tendencies (additive) 291 !! This function is called by the tendency diagnostics 291 !! This function is called by the tendency diagnostics 292 292 !! module 293 293 !!---------------------------------------------------------------------- 294 294 INTEGER , INTENT(in ) :: kt ! time step 295 295 #endif 296 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptrdx ! Temperature 296 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptrdx ! Temperature 297 297 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptrdy ! Salinity 298 298 INTEGER , INTENT(in ) :: ktrd ! tracer trend index … … 354 354 !! 355 355 !! ** Purpose : Collect momentum tendencies (additive) 356 !! This function is called by the tendency diagnostics 356 !! This function is called by the tendency diagnostics 357 357 !! module 358 358 !!---------------------------------------------------------------------- 359 359 INTEGER , INTENT(in ) :: kt ! time step 360 360 #endif 361 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptrdx ! Temperature 361 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptrdx ! Temperature 362 362 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptrdy ! Salinity 363 363 INTEGER , INTENT(in ) :: ktrd ! tracer trend index … … 469 469 !! 470 470 !! ** Purpose : Apply collinear perturbation to ice fields 471 !! For specific processes coded in LIM2/LIM3 471 !! For specific processes coded in LIM2/LIM3 472 472 !!---------------------------------------------------------------------- 473 473 ! … … 529 529 zicewrk(:,:,jm) = z5 ; jm=jm+1 530 530 zicewrk(:,:,jm) = z6 ; jm=jm+1 531 zicewrk(:,:,jm) = z7 531 zicewrk(:,:,jm) = z7 532 532 ENDIF 533 533 IF( kopt .EQ. 3 ) THEN … … 601 601 CALL lbc_lnk( u_ice, 'U', -1. ) 602 602 CALL lbc_lnk( v_ice, 'V', -1. ) 603 #endif 603 #endif 604 604 #if defined key_lim2 && defined key_lim2_vp 605 605 CALL lbc_lnk( u_ice(:,1:jpj), 'I', -1. ) 606 606 CALL lbc_lnk( v_ice(:,1:jpj), 'I', -1. ) 607 #endif 607 #endif 608 608 ENDIF 609 609 DEALLOCATE ( zicewrk ) … … 616 616 !! *** ROUTINE spp_gen *** 617 617 !! 618 !! ** Purpose : Perturbing parameters (generic function) 618 !! ** Purpose : Perturbing parameters (generic function) 619 619 !! Given a value of standard deviation, the 2D parameter 620 620 !! coeff is perturbed following an additive Normal, … … 633 633 CHARACTER (LEN=99) :: cstrng 634 634 INTEGER :: jklev 635 636 #if defined key_traldf_c3d || key_traldf_c2d 635 637 636 638 CALL wrk_alloc(jpi,jpj,gauss) … … 687 689 jklev = klev 688 690 ELSE 689 jklev = 0 691 jklev = 0 690 692 ENDIF 691 693 CALL spp_stats(kt,kspp,jklev,coeff) … … 693 695 694 696 CALL wrk_dealloc(jpi,jpj,gauss) 697 698 #else 699 CALL ctl_stop('key_traldf_c1d is not a valid key for STO') 700 #endif 695 701 696 702 END SUBROUTINE … … 707 713 REAL(wp) :: mi,ma 708 714 CHARACTER(LEN=16) :: cstr = ' ' 709 SELECT CASE ( kp ) 710 CASE( jk_spp_alb ) 715 SELECT CASE ( kp ) 716 CASE( jk_spp_alb ) 711 717 cstr = 'ALBEDO ' 712 CASE( jk_spp_rhg ) 713 cstr = 'ICE RHEOLOGY' 714 CASE( jk_spp_relw ) 715 cstr = 'RELATIVE WND' 716 CASE( jk_spp_dqdt ) 717 cstr = 'SST RELAXAT.' 718 CASE( jk_spp_deds ) 719 cstr = 'SSS RELAXAT.' 720 CASE( jk_spp_arnf ) 721 cstr = 'RIVER MIXING' 722 CASE( jk_spp_geot ) 723 cstr = 'GEOTHERM.FLX' 724 CASE( jk_spp_qsi0 ) 725 cstr = 'SOLAR EXTIN.' 726 CASE( jk_spp_bfr ) 727 cstr = 'BOTTOM FRICT' 728 CASE( jk_spp_aevd ) 729 cstr = 'EDDY VISCDIF' 730 CASE( jk_spp_avt ) 718 CASE( jk_spp_rhg ) 719 cstr = 'ICE RHEOLOGY' 720 CASE( jk_spp_relw ) 721 cstr = 'RELATIVE WND' 722 CASE( jk_spp_dqdt ) 723 cstr = 'SST RELAXAT.' 724 CASE( jk_spp_deds ) 725 cstr = 'SSS RELAXAT.' 726 CASE( jk_spp_arnf ) 727 cstr = 'RIVER MIXING' 728 CASE( jk_spp_geot ) 729 cstr = 'GEOTHERM.FLX' 730 CASE( jk_spp_qsi0 ) 731 cstr = 'SOLAR EXTIN.' 732 CASE( jk_spp_bfr ) 733 cstr = 'BOTTOM FRICT' 734 CASE( jk_spp_aevd ) 735 cstr = 'EDDY VISCDIF' 736 CASE( jk_spp_avt ) 731 737 cstr = 'VERT. DIFFUS' 732 CASE( jk_spp_avm ) 738 CASE( jk_spp_avm ) 733 739 cstr = 'VERT. VISCOS' 734 CASE( jk_spp_tkelc ) 740 CASE( jk_spp_tkelc ) 735 741 cstr = 'TKE LANGMUIR' 736 CASE( jk_spp_tkedf ) 737 cstr = 'TKE RN_EDIFF' 738 CASE( jk_spp_tkeds ) 739 cstr = 'TKE RN_EDISS' 740 CASE( jk_spp_tkebb ) 742 CASE( jk_spp_tkedf ) 743 cstr = 'TKE RN_EDIFF' 744 CASE( jk_spp_tkeds ) 745 cstr = 'TKE RN_EDISS' 746 CASE( jk_spp_tkebb ) 741 747 cstr = 'TKE RN_EBB ' 742 CASE( jk_spp_tkefr ) 748 CASE( jk_spp_tkefr ) 743 749 cstr = 'TKE RN_EFR ' 744 CASE( jk_spp_ahtu ) 750 CASE( jk_spp_ahtu ) 745 751 cstr = 'TRALDF AHTU ' 746 CASE( jk_spp_ahtv ) 752 CASE( jk_spp_ahtv ) 747 753 cstr = 'TRALDF AHTV ' 748 CASE( jk_spp_ahtw ) 754 CASE( jk_spp_ahtw ) 749 755 cstr = 'TRALDF AHTW ' 750 CASE( jk_spp_ahtt ) 756 CASE( jk_spp_ahtt ) 751 757 cstr = 'TRALDF AHTT ' 752 758 CASE( jk_spp_ahubbl ) … … 795 801 796 802 DO jp=1,jk_spp 797 SELECT CASE ( jp ) 798 CASE( jk_spp_alb ) 803 SELECT CASE ( jp ) 804 CASE( jk_spp_alb ) 799 805 cstr = 'ALBEDO ' 800 CASE( jk_spp_rhg ) 801 cstr = 'ICE RHEOLOGY' 802 CASE( jk_spp_relw ) 803 cstr = 'RELATIVE WND' 804 CASE( jk_spp_dqdt ) 805 cstr = 'SST RELAXAT.' 806 CASE( jk_spp_deds ) 807 cstr = 'SSS RELAXAT.' 808 CASE( jk_spp_arnf ) 809 cstr = 'RIVER MIXING' 810 CASE( jk_spp_geot ) 811 cstr = 'GEOTHERM.FLX' 812 CASE( jk_spp_qsi0 ) 813 cstr = 'SOLAR EXTIN.' 814 CASE( jk_spp_bfr ) 815 cstr = 'BOTTOM FRICT' 816 CASE( jk_spp_aevd ) 817 cstr = 'EDDY VISCDIF' 818 CASE( jk_spp_avt ) 806 CASE( jk_spp_rhg ) 807 cstr = 'ICE RHEOLOGY' 808 CASE( jk_spp_relw ) 809 cstr = 'RELATIVE WND' 810 CASE( jk_spp_dqdt ) 811 cstr = 'SST RELAXAT.' 812 CASE( jk_spp_deds ) 813 cstr = 'SSS RELAXAT.' 814 CASE( jk_spp_arnf ) 815 cstr = 'RIVER MIXING' 816 CASE( jk_spp_geot ) 817 cstr = 'GEOTHERM.FLX' 818 CASE( jk_spp_qsi0 ) 819 cstr = 'SOLAR EXTIN.' 820 CASE( jk_spp_bfr ) 821 cstr = 'BOTTOM FRICT' 822 CASE( jk_spp_aevd ) 823 cstr = 'EDDY VISCDIF' 824 CASE( jk_spp_avt ) 819 825 cstr = 'VERT. DIFFUS' 820 CASE( jk_spp_avm ) 826 CASE( jk_spp_avm ) 821 827 cstr = 'VERT. VISCOS' 822 CASE( jk_spp_tkelc ) 828 CASE( jk_spp_tkelc ) 823 829 cstr = 'TKE LANGMUIR' 824 CASE( jk_spp_tkedf ) 825 cstr = 'TKE RN_EDIFF' 826 CASE( jk_spp_tkeds ) 827 cstr = 'TKE RN_EDISS' 828 CASE( jk_spp_tkebb ) 830 CASE( jk_spp_tkedf ) 831 cstr = 'TKE RN_EDIFF' 832 CASE( jk_spp_tkeds ) 833 cstr = 'TKE RN_EDISS' 834 CASE( jk_spp_tkebb ) 829 835 cstr = 'TKE RN_EBB ' 830 CASE( jk_spp_tkefr ) 836 CASE( jk_spp_tkefr ) 831 837 cstr = 'TKE RN_EFR ' 832 CASE( jk_spp_ahtu ) 838 CASE( jk_spp_ahtu ) 833 839 cstr = 'TRALDF AHTU ' 834 CASE( jk_spp_ahtv ) 840 CASE( jk_spp_ahtv ) 835 841 cstr = 'TRALDF AHTV ' 836 CASE( jk_spp_ahtw ) 842 CASE( jk_spp_ahtw ) 837 843 cstr = 'TRALDF AHTW ' 838 CASE( jk_spp_ahtt ) 844 CASE( jk_spp_ahtt ) 839 845 cstr = 'TRALDF AHTT ' 840 846 CASE( jk_spp_ahubbl ) … … 922 928 INTEGER :: jk 923 929 930 #if defined key_traldf_c3d || key_traldf_c2d 931 924 932 CALL wrk_alloc(jpi,jpj,gauss) 925 933 … … 969 977 ENDIF 970 978 979 #else 980 CALL ctl_stop('key_traldf_c1d is not a valid key for STO') 981 982 #endif 983 971 984 CALL wrk_dealloc(jpi,jpj,gauss) 972 985 … … 997 1010 REAL(wp) :: zsd,xme 998 1011 INTEGER :: jk 1012 1013 #if defined key_dynldf_c3d || key_dynldf_c2d 999 1014 1000 1015 CALL wrk_alloc(jpi,jpj,gauss) … … 1046 1061 1047 1062 CALL wrk_dealloc(jpi,jpj,gauss) 1063 1064 #else 1065 CALL ctl_stop('key_traldf_c1d is not a valid key for STO') 1066 #endif 1048 1067 1049 1068 END SUBROUTINE … … 1195 1214 ! Note: sshn should be staggered before being used. 1196 1215 SELECT CASE ( cd_type ) 1197 CASE ( 'T' ) 1216 CASE ( 'T' ) 1198 1217 jk=1 1199 1218 zv = SUM( tmask_i(:,:)*tmask(:,:,jk)*e1t(:,:)*e2t(:,:)*sshn(:,:)*zts(:,:,1) ) … … 1285 1304 ! Random noise on 2d field 1286 1305 IF ( istep == 1 ) THEN 1287 CALL sppt_rand2d( g2d ) 1306 CALL sppt_rand2d( g2d ) 1288 1307 CALL lbc_lnk( g2d , 'T', 1._wp) 1289 1308 g2d_save = g2d … … 1297 1316 g2d = rn_skeb_stdev * g2d_save / rn_sppt_stdev 1298 1317 ENDIF 1299 1318 1300 1319 ! Laplacian filter and re-normalization 1301 1320 DO jf = 1, nk … … 1314 1333 ENDIF 1315 1334 #endif 1316 1335 1317 1336 ! AR(1) process and array swap 1318 1337 g2d = a*gb + b*g2d … … 1360 1379 ENDDO 1361 1380 ENDIF 1362 1381 1363 1382 ! Bound 1364 1383 IF( nn_sppt_step_bound .EQ. 2 ) THEN … … 1482 1501 1483 1502 #ifdef NEMO_V34 1484 REWIND( numnam ) 1503 REWIND( numnam ) 1485 1504 READ ( numnam, namstopack ) 1486 1505 #else 1487 REWIND( numnam_ref ) 1506 REWIND( numnam_ref ) 1488 1507 READ ( numnam_ref, namstopack, IOSTAT = ios, ERR = 901) 1489 1508 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in reference namelist', lwp ) 1490 1509 1491 REWIND( numnam_cfg ) 1510 REWIND( numnam_cfg ) 1492 1511 READ ( numnam_cfg, namstopack, IOSTAT = ios, ERR = 902 ) 1493 1512 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in configuration namelist', lwp ) … … 1568 1587 WRITE(numout,*) 1569 1588 WRITE(numout,*) ' Number of passes for spatial filter (AR1 field) spp_filter_pass:', spp_filter_pass 1570 WRITE(numout,*) ' Standard deviation of random generator (AR1 field) rn_spp_stdev :', rn_spp_stdev 1589 WRITE(numout,*) ' Standard deviation of random generator (AR1 field) rn_spp_stdev :', rn_spp_stdev 1571 1590 WRITE(numout,*) ' Decorr. time scale (AR1 field) rn_spp_tau :', rn_spp_tau 1572 1591 WRITE(numout,*) 1573 WRITE(numout,*) ' SPP for bottom friction coeff nn_spp_bfr :', nn_spp_bfr 1574 WRITE(numout,*) ' STDEV rn_bfr_sd :', rn_bfr_sd 1575 WRITE(numout,*) ' SPP for SST relaxation coeff nn_spp_dqdt :', nn_spp_dqdt 1576 WRITE(numout,*) ' STDEV rn_dqdt_sd :', rn_dqdt_sd 1577 WRITE(numout,*) ' SPP for SSS relaxation coeff nn_spp_dedt :', nn_spp_dedt 1578 WRITE(numout,*) ' STDEV rn_dedt_sd :', rn_dedt_sd 1579 WRITE(numout,*) ' SPP for vertical tra mixing coeff (only TKE, GLS) nn_spp_avt :', nn_spp_avt 1580 WRITE(numout,*) ' STDEV rn_avt_sd :', rn_avt_sd 1581 WRITE(numout,*) ' SPP for vertical dyn mixing coeff (only TKE, GLS) nn_spp_avm :', nn_spp_avm 1582 WRITE(numout,*) ' STDEV rn_avm_sd :', rn_avm_sd 1592 WRITE(numout,*) ' SPP for bottom friction coeff nn_spp_bfr :', nn_spp_bfr 1593 WRITE(numout,*) ' STDEV rn_bfr_sd :', rn_bfr_sd 1594 WRITE(numout,*) ' SPP for SST relaxation coeff nn_spp_dqdt :', nn_spp_dqdt 1595 WRITE(numout,*) ' STDEV rn_dqdt_sd :', rn_dqdt_sd 1596 WRITE(numout,*) ' SPP for SSS relaxation coeff nn_spp_dedt :', nn_spp_dedt 1597 WRITE(numout,*) ' STDEV rn_dedt_sd :', rn_dedt_sd 1598 WRITE(numout,*) ' SPP for vertical tra mixing coeff (only TKE, GLS) nn_spp_avt :', nn_spp_avt 1599 WRITE(numout,*) ' STDEV rn_avt_sd :', rn_avt_sd 1600 WRITE(numout,*) ' SPP for vertical dyn mixing coeff (only TKE, GLS) nn_spp_avm :', nn_spp_avm 1601 WRITE(numout,*) ' STDEV rn_avm_sd :', rn_avm_sd 1583 1602 WRITE(numout,*) ' SPP for solar penetration scheme (only RGB) nn_spp_qsi0 :', nn_spp_qsi0 1584 WRITE(numout,*) ' STDEV rn_qsi0_sd :', rn_qsi0_sd 1603 WRITE(numout,*) ' STDEV rn_qsi0_sd :', rn_qsi0_sd 1585 1604 WRITE(numout,*) ' SPP for horiz. diffusivity U nn_spp_ahtu :', nn_spp_ahtu 1586 WRITE(numout,*) ' STDEV rn_ahtu_sd :', rn_ahtu_sd 1605 WRITE(numout,*) ' STDEV rn_ahtu_sd :', rn_ahtu_sd 1587 1606 WRITE(numout,*) ' SPP for horiz. diffusivity V nn_spp_ahtv :', nn_spp_ahtv 1588 WRITE(numout,*) ' STDEV rn_ahtv_sd :', rn_ahtv_sd 1607 WRITE(numout,*) ' STDEV rn_ahtv_sd :', rn_ahtv_sd 1589 1608 WRITE(numout,*) ' SPP for horiz. diffusivity W nn_spp_ahtw :', nn_spp_ahtw 1590 WRITE(numout,*) ' STDEV rn_ahtw_sd :', rn_ahtw_sd 1609 WRITE(numout,*) ' STDEV rn_ahtw_sd :', rn_ahtw_sd 1591 1610 WRITE(numout,*) ' SPP for horiz. diffusivity T nn_spp_ahtt :', nn_spp_ahtt 1592 WRITE(numout,*) ' STDEV rn_ahtt_sd :', rn_ahtt_sd 1611 WRITE(numout,*) ' STDEV rn_ahtt_sd :', rn_ahtt_sd 1593 1612 WRITE(numout,*) ' SPP for horiz. viscosity (1/3) nn_spp_ahm1 :', nn_spp_ahm1 1594 WRITE(numout,*) ' STDEV rn_ahm1_sd :', rn_ahm1_sd 1613 WRITE(numout,*) ' STDEV rn_ahm1_sd :', rn_ahm1_sd 1595 1614 WRITE(numout,*) ' SPP for horiz. viscosity (2/4) nn_spp_ahm2 :', nn_spp_ahm2 1596 WRITE(numout,*) ' STDEV rn_ahm2_sd :', rn_ahm2_sd 1615 WRITE(numout,*) ' STDEV rn_ahm2_sd :', rn_ahm2_sd 1597 1616 WRITE(numout,*) ' SPP for relative wind factor nn_spp_relw :', nn_spp_relw 1598 1617 WRITE(numout,*) ' (use 4, 5, 6 for nn_spp_relw to have options 1, 2, 3 with limits bounded to [0,1]' 1599 WRITE(numout,*) ' STDEV rn_relw_sd :', rn_relw_sd 1618 WRITE(numout,*) ' STDEV rn_relw_sd :', rn_relw_sd 1600 1619 WRITE(numout,*) ' SPP for mixing close to river mouth nn_spp_arnf :', nn_spp_arnf 1601 WRITE(numout,*) ' STDEV rn_arnf_sd :', rn_arnf_sd 1620 WRITE(numout,*) ' STDEV rn_arnf_sd :', rn_arnf_sd 1602 1621 WRITE(numout,*) ' SPP for geothermal heating nn_spp_geot :', nn_spp_geot 1603 WRITE(numout,*) ' STDEV rn_geot_sd :', rn_geot_sd 1622 WRITE(numout,*) ' STDEV rn_geot_sd :', rn_geot_sd 1604 1623 WRITE(numout,*) ' SPP for enhanced vertical diffusion nn_spp_aevd :', nn_spp_aevd 1605 WRITE(numout,*) ' STDEV rn_aevd_sd :', rn_aevd_sd 1624 WRITE(numout,*) ' STDEV rn_aevd_sd :', rn_aevd_sd 1606 1625 WRITE(numout,*) ' SPP for TKE rn_lc Langmuir cell coefficient nn_spp_tkelc :', nn_spp_tkelc 1607 WRITE(numout,*) ' STDEV rn_tkelc_sd :', rn_tkelc_sd 1626 WRITE(numout,*) ' STDEV rn_tkelc_sd :', rn_tkelc_sd 1608 1627 WRITE(numout,*) ' SPP for TKE rn_ediff Eddy diff. coefficient nn_spp_tkedf :', nn_spp_tkedf 1609 WRITE(numout,*) ' STDEV rn_tkedf_sd :', rn_tkedf_sd 1628 WRITE(numout,*) ' STDEV rn_tkedf_sd :', rn_tkedf_sd 1610 1629 WRITE(numout,*) ' SPP for TKE rn_ediss Kolmogoroff dissipation coeff. nn_spp_tkeds :', nn_spp_tkeds 1611 WRITE(numout,*) ' STDEV rn_tkeds_sd :', rn_tkeds_sd 1630 WRITE(numout,*) ' STDEV rn_tkeds_sd :', rn_tkeds_sd 1612 1631 WRITE(numout,*) ' SPP for TKE rn_ebb Surface input of tke nn_spp_tkebb :', nn_spp_tkebb 1613 WRITE(numout,*) ' STDEV rn_tkebb_sd :', rn_tkebb_sd 1632 WRITE(numout,*) ' STDEV rn_tkebb_sd :', rn_tkebb_sd 1614 1633 WRITE(numout,*) ' SPP for TKE rn_efr Fraction of srf TKE below ML nn_spp_tkefr :', nn_spp_tkefr 1615 WRITE(numout,*) ' STDEV rn_tkefr_sd :', rn_tkefr_sd 1634 WRITE(numout,*) ' STDEV rn_tkefr_sd :', rn_tkefr_sd 1616 1635 WRITE(numout,*) ' SPP for BBL U diffusivity nn_spp_ahubbl:', nn_spp_ahubbl 1617 1636 WRITE(numout,*) ' STDEV rn_ahubbl_sd :', rn_ahubbl_sd … … 1626 1645 WRITE(numout,*) 1627 1646 WRITE(numout,*) ' SKEB Perturbation scheme ' 1628 WRITE(numout,*) ' SKEB switch ln_skeb :', ln_skeb 1629 WRITE(numout,*) ' SKEB ratio of backscattered energy rn_skeb :', rn_skeb 1647 WRITE(numout,*) ' SKEB switch ln_skeb :', ln_skeb 1648 WRITE(numout,*) ' SKEB ratio of backscattered energy rn_skeb :', rn_skeb 1630 1649 WRITE(numout,*) ' Frequency update for dissipation mask nn_dcom_freq :', nn_dcom_freq 1631 1650 WRITE(numout,*) ' Numerical dissipation factor (resolut. dependent) rn_kh :', rn_kh 1632 1651 WRITE(numout,*) ' Number of passes for spatial filter (AR1 field) skeb_filter_pass:', skeb_filter_pass 1633 WRITE(numout,*) ' Standard deviation of random generator (AR1 field) rn_skeb_stdev:', rn_skeb_stdev 1652 WRITE(numout,*) ' Standard deviation of random generator (AR1 field) rn_skeb_stdev:', rn_skeb_stdev 1634 1653 WRITE(numout,*) ' Decorr. time scale (AR1 field) rn_skeb_tau :', rn_skeb_tau 1635 1654 WRITE(numout,*) ' Option of convection energy dissipation nn_dconv :', nn_dconv … … 1752 1771 1753 1772 ! Find filter attenuation factor 1754 1773 1755 1774 flt_fac = sppt_fltfac( sppt_filter_pass ) 1756 1775 rdt_sppt = nn_rndm_freq * rn_rdt 1757 1776 1758 1777 IF( ln_sppt_taumap ) THEN 1759 1778 CALL iom_open ( 'sppt_tau_map', inum ) … … 1798 1817 gauss_b = 0._wp 1799 1818 ! Weigths 1800 gauss_w(:) = 1.0_wp 1819 gauss_w(:) = 1.0_wp 1801 1820 IF( nn_vwei .eq. 1 ) THEN 1802 1821 gauss_w(1) = 0.0_wp … … 1861 1880 IF(lwp .and. ln_stopack_diags) & 1862 1881 CALL ctl_opn(numdiag, 'stopack.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 1863 1882 1864 1883 END SUBROUTINE stopack_init 1865 1884 ! … … 1874 1893 INTEGER :: id1, jseed 1875 1894 CHARACTER(LEN=10) :: clseed='spsd0_0000' 1876 INTEGER(KIND=8) :: ziseed(4) ! RNG seeds in integer type 1877 INTEGER(KIND=8) :: ivals(8) 1895 INTEGER(KIND=8) :: ziseed(4) ! RNG seeds in integer type 1896 INTEGER(KIND=8) :: ivals(8) 1878 1897 REAL(wp) :: zrseed4(4) ! RNG seeds in integer type 1879 1898 REAL(wp) :: zrseed2d(jpi,jpj) … … 1983 2002 !!--------------------------------------------------------------------- 1984 2003 ! 1985 ALLOCATE( spptt(jpi,jpj,jpk) , sppts(jpi,jpj,jpk) , gauss_n(jpi,jpj,jpk) ,& 1986 gauss_nu(jpi,jpj,jpk) , gauss_nv(jpi,jpj,jpk) , & 2004 ALLOCATE( spptt(jpi,jpj,jpk) , sppts(jpi,jpj,jpk) , gauss_n(jpi,jpj,jpk) ,& 2005 gauss_nu(jpi,jpj,jpk) , gauss_nv(jpi,jpj,jpk) , & 1987 2006 spptu(jpi,jpj,jpk) , spptv(jpi,jpj,jpk) , gauss_n_2d(jpi,jpj) ,& 1988 2007 gauss_b (jpi,jpj), sppt_tau(jpi,jpj), sppt_a(jpi,jpj), sppt_b(jpi,jpj), gauss_w(jpk),& … … 2208 2227 IF ( ln_dpsiv ) THEN 2209 2228 DO jp=1,jpni-1 2210 IF( jpri == jp ) THEN ! SEND TO EAST 2229 IF( jpri == jp ) THEN ! SEND TO EAST 2211 2230 zwrk(1:jpj) = dpsiv(jpi-1,:) 2212 2231 tag=2000+narea … … 2268 2287 REAL :: ds,dt,dtot,kh2 2269 2288 INTEGER :: ji,jj,jk 2270 2289 2271 2290 IF ( mt .eq. nit000 ) THEN 2272 2291 ALLOCATE ( dnum(jpi,jpj,jpk) ) … … 2287 2306 dt = (vn(ji,jj,jk)-vn(ji-1,jj,jk))*vmask(ji,jj,jk)*vmask(ji-1,jj,jk)/ e2v(ji,jj) + & 2288 2307 (un(ji,jj,jk)-un(ji,jj-1,jk))*umask(ji,jj,jk)*umask(ji,jj-1,jk)/ e1u(ji,jj) 2289 2308 2290 2309 dtot = sqrt( ds*ds + dt*dt ) * tmask(ji,jj,jk) 2291 2310 dnum(ji,jj,jk) = dtot*dtot*dtot*kh2*e1t(ji,jj)*e2t(ji,jj) … … 2293 2312 ENDDO 2294 2313 ENDDO 2295 2314 2296 2315 CALL lbc_lnk(dnum,'T',1._wp) 2297 2316 2298 2317 ENDIF 2299 2318 2300 END SUBROUTINE 2319 END SUBROUTINE 2301 2320 2302 2321 SUBROUTINE skeb_dcon ( mt ) … … 2329 2348 2330 2349 zz = - grav*avt(ji,jj,jk) * ( rhd(ji,jj,jk)-rhd(ji,jj,jk-1) ) * wmask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) & 2331 & / ( rau0 * fse3w(ji,jj,jk) ) 2350 & / ( rau0 * fse3w(ji,jj,jk) ) 2332 2351 2333 2352 dcon(ji,jj,jk) = kc2*zz*e1t(ji,jj)*e2t(ji,jj)*rau0 / fse3w(ji,jj,jk) … … 2378 2397 IF(ln_skeb_own_gauss) THEN 2379 2398 DO jk=1,jpkm1 2380 psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) * gauss_n_2d_k(:,:) 2399 psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) * gauss_n_2d_k(:,:) 2381 2400 ENDDO 2382 2401 ELSE … … 2407 2426 IF(ln_skeb_own_gauss) THEN 2408 2427 DO jk=1,jpkm1 2409 psi(:,:,jk) = rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 2428 psi(:,:,jk) = rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 2410 2429 ENDDO 2411 2430 ELSE … … 2440 2459 IF(ln_skeb_own_gauss) THEN 2441 2460 DO jk=1,jpkm1 2442 psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 2461 psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 2443 2462 ENDDO 2444 2463 ELSE -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
r11442 r13355 12 12 13 13 !!---------------------------------------------------------------------- 14 !! tra_bbc : update the tracer trend at ocean bottom 14 !! tra_bbc : update the tracer trend at ocean bottom 15 15 !! tra_bbc_init : initialization of geothermal heat flux trend 16 16 !!---------------------------------------------------------------------- … … 19 19 USE phycst ! physical constants 20 20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! trends manager: tracers 21 USE trdtra ! trends manager: tracers 22 22 USE in_out_manager ! I/O manager 23 23 USE iom ! I/O manager … … 44 44 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: qgh_trd1 ! geothermal heating trend 45 45 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_qgh ! structure of input qgh (file informations, fields read) 46 46 47 47 !! * Substitutions 48 48 # include "domzgr_substitute.h90" … … 58 58 !! *** ROUTINE tra_bbc *** 59 59 !! 60 !! ** Purpose : Compute the bottom boundary contition on temperature 61 !! associated with geothermal heating and add it to the 60 !! ** Purpose : Compute the bottom boundary contition on temperature 61 !! associated with geothermal heating and add it to the 62 62 !! general trend of temperature equations. 63 63 !! 64 !! ** Method : The geothermal heat flux set to its constant value of 64 !! ** Method : The geothermal heat flux set to its constant value of 65 65 !! 86.4 mW/m2 (Stein and Stein 1992, Huang 1999). 66 66 !! The temperature trend associated to this heat flux through the … … 92 92 ! ! Add the geothermal heat flux trend on temperature 93 93 94 #if defined key_traldf_c2d || key_traldf_c3d 94 95 IF( ln_stopack .AND. nn_spp_geot > 0) THEN 95 96 qgh_trd1(:,:) = qgh_trd0(:,:) 96 97 CALL spp_gen(kt, qgh_trd1, nn_spp_geot, rn_geot_sd, jk_spp_geot) 97 98 ENDIF 99 #else 100 IF ( ln_stopack .AND. nn_spp_geot > 0 ) & 101 & CALL ctl_stop( 'tra_bbc: parameter perturbation will only work with '// & 102 'key_traldf_c2d or key_traldf_c3d') 103 #endif 104 105 98 106 DO jj = 2, jpjm1 99 107 DO ji = 2, jpim1 … … 144 152 CHARACTER(len=256) :: cn_dir ! Root directory for location of ssr files 145 153 ! 146 NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 154 NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 147 155 !!---------------------------------------------------------------------- 148 156 -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
r11442 r13355 32 32 USE trdtra ! trends: active tracers 33 33 ! 34 USE iom ! IOM library 34 USE iom ! IOM library 35 35 USE in_out_manager ! I/O manager 36 36 USE lbclnk ! ocean lateral boundary conditions … … 38 38 USE wrk_nemo ! Memory Allocation 39 39 USE timing ! Timing 40 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 40 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 41 41 USE stopack 42 42 … … 199 199 ! 200 200 ahu_bbl_1(:,:) = ahu_bbl(:,:) 201 #if defined key_traldf_c2d || key_traldf_c3d 201 202 IF( ln_stopack .AND. nn_spp_ahubbl > 0 ) THEN 202 203 CALL spp_gen(1, ahu_bbl_1, nn_spp_ahubbl, rn_ahubbl_sd, jk_spp_ahubbl ) 203 204 ENDIF 205 #else 206 IF ( ln_stopack .AND. nn_spp_ahubbl > 0 ) & 207 & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// & 208 'key_traldf_c2d or key_traldf_c3d') 209 #endif 210 211 204 212 ahv_bbl_1(:,:) = ahv_bbl(:,:) 213 #if defined key_traldf_c2d || key_traldf_c3d 205 214 IF( ln_stopack .AND. nn_spp_ahvbbl > 0 ) THEN 206 215 CALL spp_gen(1, ahv_bbl_1, nn_spp_ahvbbl, rn_ahvbbl_sd, jk_spp_ahvbbl ) 207 216 ENDIF 217 #else 218 IF ( ln_stopack .AND. nn_spp_ahvbbl > 0 ) & 219 & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// & 220 'key_traldf_c2d or key_traldf_c3d') 221 #endif 222 208 223 ! 209 224 DO jn = 1, kjpt ! tracer loop … … 215 230 END DO 216 231 END DO 217 ! 232 ! 218 233 DO jj = 2, jpjm1 ! Compute the trend 219 234 DO ji = 2, jpim1 … … 426 441 za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point 427 442 zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 428 ! ! 2*masked bottom density gradient 443 ! ! 2*masked bottom density gradient 429 444 zgdrho = ( za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) ) & 430 445 - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) ) ) * umask(ji,jj,1) … … 585 600 gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) ) ) 586 601 ENDIF 587 ! 602 ! 588 603 IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN 589 604 mgrhv(ji,jj) = INT( SIGN( 1.e0, & -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r11442 r13355 9 9 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 10 10 !! - ! 2005-11 (G. Madec) zco, zps, sco coordinate 11 !! 3.2 ! 2009-04 (G. Madec & NEMO team) 12 !! 3.4 ! 2012-05 (C. Rousset) store attenuation coef for use in ice model 11 !! 3.2 ! 2009-04 (G. Madec & NEMO team) 12 !! 3.4 ! 2012-05 (C. Rousset) store attenuation coef for use in ice model 13 13 !! 3.6 ! 2015-12 (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll 14 14 !!---------------------------------------------------------------------- … … 43 43 ! !!* Namelist namtra_qsr: penetrative solar radiation 44 44 LOGICAL , PUBLIC :: ln_traqsr !: light absorption (qsr) flag 45 LOGICAL , PUBLIC :: ln_qsr_rgb !: Red-Green-Blue light absorption flag 45 LOGICAL , PUBLIC :: ln_qsr_rgb !: Red-Green-Blue light absorption flag 46 46 LOGICAL , PUBLIC :: ln_qsr_2bd !: 2 band light absorption flag 47 47 LOGICAL , PUBLIC :: ln_qsr_bio !: bio-model light absorption flag … … 51 51 REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands) 52 52 REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands) 53 53 54 54 ! Module variables 55 55 REAL(wp), ALLOCATABLE :: xsi0r(:,:) !: inverse of rn_si0 … … 80 80 !! Considering the 2 wavebands case: 81 81 !! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 82 !! The temperature trend associated with the solar radiation penetration 82 !! The temperature trend associated with the solar radiation penetration 83 83 !! is given by : zta = 1/e3t dk[ I ] / (rau0*Cp) 84 84 !! At the bottom, boudary condition for the radiation is no flux : … … 86 86 !! in the last ocean level. 87 87 !! In z-coordinate case, the computation is only done down to the 88 !! level where I(k) < 1.e-15 W/m2. In addition, the coefficients 88 !! level where I(k) < 1.e-15 W/m2. In addition, the coefficients 89 89 !! used for the computation are calculated one for once as they 90 90 !! depends on k only. … … 106 106 REAL(wp) :: zz0, zz1, z1_e3t ! - - 107 107 REAL(wp) :: zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze 108 REAL(wp) :: zlogc, zlogc2, zlogc3 108 REAL(wp) :: zlogc, zlogc2, zlogc3 109 109 REAL(wp), POINTER, DIMENSION(:,: ) :: zekb, zekg, zekr 110 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt, zchl3d … … 113 113 IF( nn_timing == 1 ) CALL timing_start('tra_qsr') 114 114 ! 115 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 116 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 115 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 116 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 117 117 ! 118 118 IF( kt == nit000 ) THEN … … 124 124 125 125 IF( l_trdtra ) THEN ! Save ta and sa trends 126 CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 126 CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 127 127 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 128 128 ENDIF … … 150 150 ! Compute now qsr tracer content field 151 151 ! ************************************ 152 152 153 153 ! ! ============================================== ! 154 154 IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates ! … … 159 159 ! Add to the general trend 160 160 DO jk = 1, jpkm1 161 DO jj = 2, jpjm1 161 DO jj = 2, jpjm1 162 162 DO ji = fs_2, fs_jpim1 ! vector opt. 163 163 z1_e3t = zfact / fse3t(ji,jj,jk) … … 180 180 ENDIF 181 181 ! ! ============================================== ! 182 ELSE ! Ocean alone : 182 ELSE ! Ocean alone : 183 183 ! ! ============================================== ! 184 184 ! 185 ! 185 ! 186 #if defined key_traldf_c2d || key_traldf_c3d 186 187 IF( ln_stopack .AND. ( nn_spp_qsi0 > 0 ) ) THEN 187 xsi0r = rn_si0 188 CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 ) 189 xsi0r = 1.e0 / xsi0r 190 ENDIF 188 xsi0r = rn_si0 189 CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 ) 190 xsi0r = 1.e0 / xsi0r 191 ENDIF 192 #else 193 IF ( ln_stopack .AND. nn_spp_qsi0 > 0 ) & 194 & CALL ctl_stop( 'tra_qsr: parameter perturbation will only work with '// & 195 'key_traldf_c2d or key_traldf_c3d') 196 #endif 197 191 198 ! ! ------------------------- ! 192 199 IF( ln_qsr_rgb) THEN ! R-G-B light penetration ! … … 199 206 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 200 207 DO jk = 1, nksr + 1 201 zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 208 zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 202 209 ENDDO 203 210 ! … … 220 227 zpsimax = 0.6 - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3 221 228 zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2 222 zCze = 1.12 * (zchl)**0.803 229 zCze = 1.12 * (zchl)**0.803 223 230 DO jk = 1, nksr + 1 224 231 zpsi = fsdept(ji,jj,jk) / zze … … 231 238 ELSE !* Variable ocean volume but constant chrlorophyll 232 239 DO jk = 1, nksr + 1 233 zchl3d(:,:,jk) = 0.05 240 zchl3d(:,:,jk) = 0.05 234 241 ENDDO 235 242 ENDIF … … 256 263 !CDIR NOVERRCHK 257 264 DO jj = 1, jpj 258 !CDIR NOVERRCHK 265 !CDIR NOVERRCHK 259 266 DO ji = 1, jpi 260 267 zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r(ji,jj) ) … … 289 296 END DO 290 297 END DO 291 ! 298 ! 292 299 DO jj = 1, jpj 293 300 DO ji = 1, jpi … … 296 303 zc2 = zcoef * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 297 304 zc3 = zcoef * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 298 fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,2) 305 fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,2) 299 306 END DO 300 307 END DO … … 321 328 zz0 = rn_abs * r1_rau0_rcp 322 329 zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 323 DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m 330 DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m 324 331 DO jj = 1, jpj 325 332 DO ji = 1, jpi 326 333 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r ) 327 334 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 328 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 335 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 329 336 END DO 330 337 END DO … … 344 351 DO jj = 2, jpjm1 345 352 DO ji = fs_2, fs_jpim1 ! vector opt. 346 ! (ISF) no light penetration below the ice shelves 353 ! (ISF) no light penetration below the ice shelves 347 354 qsr_hc(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1) 348 355 END DO … … 360 367 ! Add to the general trend 361 368 DO jk = 1, nksr 362 DO jj = 2, jpjm1 369 DO jj = 2, jpjm1 363 370 DO ji = fs_2, fs_jpim1 ! vector opt. 364 371 z1_e3t = zfact / fse3t(ji,jj,jk) … … 376 383 & 'at it= ', kt,' date= ', ndastp 377 384 IF(lwp) WRITE(numout,*) '~~~~' 378 IF(nn_timing == 2) CALL timing_start('iom_rstput') 385 IF(nn_timing == 2) CALL timing_start('iom_rstput') 379 386 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b' , qsr_hc ) 380 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) ! default definition in sbcssm 387 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) ! default definition in sbcssm 381 388 IF(nn_timing == 2) CALL timing_stop('iom_rstput') 382 389 ! … … 386 393 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 387 394 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 388 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 395 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 389 396 ENDIF 390 397 ! ! print mean trends (used for debugging) 391 398 IF(ln_ctl) CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 392 399 ! 393 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 394 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 400 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 401 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 395 402 ! 396 403 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr') … … 408 415 !! from two length scale of penetration (rn_si0,rn_si1) and a ratio 409 416 !! (rn_abs). These parameters are read in the namtra_qsr namelist. The 410 !! default values correspond to clear water (type I in Jerlov' 417 !! default values correspond to clear water (type I in Jerlov' 411 418 !! (1968) classification. 412 419 !! called by tra_qsr at the first timestep (nit000) … … 435 442 IF( nn_timing == 1 ) CALL timing_start('tra_qsr_init') 436 443 ! 437 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 438 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 444 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 445 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 439 446 ! 440 447 … … 465 472 466 473 IF( ln_traqsr ) THEN ! control consistency 467 ! 474 ! 468 475 IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio ) THEN 469 476 CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' ) … … 480 487 & ' 2 bands, 3 RGB bands or bio-model light penetration' ) 481 488 ! 482 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1 489 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1 483 490 IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = 2 484 491 IF( ln_qsr_rgb .AND. nn_chldta == 2 ) nqsr = 3 … … 497 504 ENDIF 498 505 ! ! ===================================== ! 499 IF( ln_traqsr ) THEN ! Initialisation of Light Penetration ! 506 IF( ln_traqsr ) THEN ! Initialisation of Light Penetration ! 500 507 ! ! ===================================== ! 501 508 ! … … 539 546 zchl = 0.05 ! constant chlorophyll 540 547 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 541 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 548 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 542 549 zekg(:,:) = rkrgb(2,irgb) 543 550 zekr(:,:) = rkrgb(3,irgb) … … 546 553 ze0(:,:,1) = rn_abs 547 554 ze1(:,:,1) = zcoef 548 ze2(:,:,1) = zcoef 555 ze2(:,:,1) = zcoef 549 556 ze3(:,:,1) = zcoef 550 557 zea(:,:,1) = tmask(:,:,1) ! = ( ze0+ze1+z2+ze3 ) * tmask 551 558 552 559 DO jk = 2, nksr+1 553 560 !CDIR NOVERRCHK 554 561 DO jj = 1, jpj 555 !CDIR NOVERRCHK 562 !CDIR NOVERRCHK 556 563 DO ji = 1, jpi 557 564 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r(ji,jj) ) … … 566 573 END DO 567 574 END DO 568 END DO 575 END DO 569 576 ! 570 577 DO jk = 1, nksr … … 598 605 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r ) 599 606 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 600 etot3(ji,jj,jk) = ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) * tmask(ji,jj,1) 607 etot3(ji,jj,jk) = ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) * tmask(ji,jj,1) 601 608 END DO 602 609 END DO … … 607 614 ENDIF 608 615 ! ! ===================================== ! 609 ELSE ! No light penetration ! 616 ELSE ! No light penetration ! 610 617 ! ! ===================================== ! 611 618 IF(lwp) THEN … … 625 632 ENDIF 626 633 ! 627 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 628 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 634 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 635 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 629 636 ! 630 637 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_init') -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r11442 r13355 107 107 WRITE(numout,*) '~~~~~~~~' 108 108 ENDIF 109 ! 109 ! 110 #if defined key_traldf_c2d || key_traldf_c3d 110 111 IF( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) THEN 111 112 bfrcoef2d = bfrcoef2d0 112 113 CALL spp_gen(kt, bfrcoef2d, nn_spp_bfr, rn_bfr_sd, jk_spp_bfr) 113 114 ENDIF 115 #else 116 IF ( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) & 117 & CALL ctl_stop( 'zdf_bfr: parameter perturbation will only work with '// & 118 'key_traldf_c2d or key_traldf_c3d') 119 #endif 114 120 ! 115 121 IF( nn_bfr == 2 ) THEN ! quadratic bottom friction only … … 140 146 END DO 141 147 END IF 142 ! 148 ! 143 149 ELSE 144 150 zbfrt(:,:) = bfrcoef2d(:,:) … … 183 189 DO ji = 2, jpim1 184 190 ! (ISF) ======================================================================== 185 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 191 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 186 192 ikbv = mikv(ji,jj) ! (1st wet ocean u- and v-points) 187 193 ! … … 363 369 bfrcoef2d(:,:) = rn_bfri2 ! initialize bfrcoef2d to the namelist variable 364 370 ENDIF 365 371 366 372 IF ( ln_isfcav ) THEN 367 373 IF(ln_tfr2d) THEN -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
r11442 r13355 20 20 USE zdfkpp ! KPP vertical mixing 21 21 USE trd_oce ! trends: ocean variables 22 USE trdtra ! trends manager: tracers 22 USE trdtra ! trends manager: tracers 23 23 USE in_out_manager ! I/O manager 24 24 USE iom ! for iom_put 25 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 26 26 USE timing ! Timing 27 USE stopack 27 USE stopack 28 28 29 29 IMPLICIT NONE … … 45 45 !!---------------------------------------------------------------------- 46 46 !! *** ROUTINE zdf_evd *** 47 !! 47 !! 48 48 !! ** Purpose : Local increased the vertical eddy viscosity and diffu- 49 49 !! sivity coefficients when a static instability is encountered. 50 50 !! 51 51 !! ** Method : avt, avm, and the 4 neighbouring avmu, avmv coefficients 52 !! are set to avevd (namelist parameter) if the water column is 52 !! are set to avevd (namelist parameter) if the water column is 53 53 !! statically unstable (i.e. if rn2 < -1.e-12 ) 54 54 !! … … 77 77 zavt_evd(:,:,:) = avt(:,:,:) ! set avt prior to evd application 78 78 79 #if defined key_traldf_c2d || key_traldf_c3d 79 80 IF( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) THEN 80 81 rn_avevd0(:,:) = rn_avevd 81 82 CALL spp_gen(kt, rn_avevd0, nn_spp_aevd, rn_aevd_sd, jk_spp_aevd) 82 83 ENDIF 84 #else 85 IF ( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) & 86 & CALL ctl_stop( 'zdf_evd: parameter perturbation will only work with '// & 87 'key_traldf_c2d or key_traldf_c3d') 88 #endif 83 89 84 90 SELECT CASE ( nn_evdm ) … … 88 94 zavm_evd(:,:,:) = avm(:,:,:) ! set avm prior to evd application 89 95 ! 90 DO jk = 1, jpkm1 96 DO jk = 1, jpkm1 91 97 DO jj = 2, jpj ! no vector opt. 92 98 DO ji = 2, jpi … … 106 112 END DO 107 113 END DO 108 END DO 114 END DO 109 115 CALL lbc_lnk( avt , 'W', 1. ) ; CALL lbc_lnk( avm , 'W', 1. ) ! Lateral boundary conditions 110 116 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) … … 113 119 CALL iom_put( "avm_evd", zavm_evd ) ! output this change 114 120 ! 115 CASE DEFAULT ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 121 CASE DEFAULT ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 116 122 DO jk = 1, jpkm1 117 !!! WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd ! agissant sur T SEUL! 123 !!! WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd ! agissant sur T SEUL! 118 124 DO jj = 1, jpj ! loop over the whole domain (no lbc_lnk call) 119 125 DO ji = 1, jpi 120 126 #if defined key_zdfkpp 121 127 ! no evd mixing in the boundary layer with KPP 122 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 .AND. fsdepw(ji,jj,jk) > hkpp(ji,jj) ) & 128 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 .AND. fsdepw(ji,jj,jk) > hkpp(ji,jj) ) & 123 129 #else 124 130 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) & … … 129 135 END DO 130 136 ! 131 END SELECT 137 END SELECT 132 138 133 139 zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:) ! change in avt due to evd -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r11442 r13355 2 2 !!====================================================================== 3 3 !! *** MODULE zdfgls *** 4 !! Ocean physics: vertical mixing coefficient computed from the gls 4 !! Ocean physics: vertical mixing coefficient computed from the gls 5 5 !! turbulent closure parameterization 6 6 !!====================================================================== … … 16 16 !! gls_rst : read/write gls restart in ocean restart file 17 17 !!---------------------------------------------------------------------- 18 USE oce ! ocean dynamics and active tracers 18 USE oce ! ocean dynamics and active tracers 19 19 USE dom_oce ! ocean space and time domain 20 20 USE domvvl ! ocean space and time domain : variable volume layer … … 31 31 USE iom ! I/O manager library 32 32 USE timing ! Timing 33 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 33 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 34 34 USE stopack 35 35 … … 62 62 REAL(wp) :: rn_crban ! Craig and Banner constant for surface breaking waves mixing 63 63 REAL(wp) :: rn_hsro ! Minimum surface roughness 64 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 64 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 65 65 66 66 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters 67 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 68 REAL(wp) :: rl_sf = 0.2_wp ! 0 <rl_sf<vkarmn 67 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 68 REAL(wp) :: rl_sf = 0.2_wp ! 0 <rl_sf<vkarmn 69 69 REAL(wp) :: rghmin = -0.28_wp 70 70 REAL(wp) :: rgh0 = 0.0329_wp … … 73 73 REAL(wp) :: ra2 = 0.74_wp 74 74 REAL(wp) :: rb1 = 16.60_wp 75 REAL(wp) :: rb2 = 10.10_wp 76 REAL(wp) :: re2 = 1.33_wp 75 REAL(wp) :: rb2 = 10.10_wp 76 REAL(wp) :: re2 = 1.33_wp 77 77 REAL(wp) :: rl1 = 0.107_wp 78 78 REAL(wp) :: rl2 = 0.0032_wp … … 134 134 INTEGER :: ji, jj, jk, ibot, ibotm1, dir ! dummy loop arguments 135 135 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1, zex2 ! local scalars 136 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 136 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 137 137 REAL(wp) :: zratio, zrn2, zflxb, sh ! - - 138 138 REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - … … 140 140 REAL(wp), POINTER, DIMENSION(:,: ) :: zdep 141 141 REAL(wp), POINTER, DIMENSION(:,: ) :: zkar 142 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 142 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 143 143 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) 144 144 REAL(wp), POINTER, DIMENSION(:,:,:) :: eb ! tke at time before … … 157 157 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 158 158 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 159 159 160 160 ! Preliminary computing 161 161 … … 166 166 avm (:,:,:) = avm_k (:,:,:) 167 167 avmu(:,:,:) = avmu_k(:,:,:) 168 avmv(:,:,:) = avmv_k(:,:,:) 168 avmv(:,:,:) = avmv_k(:,:,:) 169 169 ENDIF 170 170 171 171 ! Compute surface and bottom friction at T-points 172 !CDIR NOVERRCHK 173 DO jj = 2, jpjm1 174 !CDIR NOVERRCHK 175 DO ji = fs_2, fs_jpim1 ! vector opt. 172 !CDIR NOVERRCHK 173 DO jj = 2, jpjm1 174 !CDIR NOVERRCHK 175 DO ji = fs_2, fs_jpim1 ! vector opt. 176 176 ! 177 177 ! surface friction 178 178 ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 179 ! 180 ! bottom friction (explicit before friction) 181 ! Note that we chose here not to bound the friction as in dynbfr) 182 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 183 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 184 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 185 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 186 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 187 END DO 188 END DO 179 ! 180 ! bottom friction (explicit before friction) 181 ! Note that we chose here not to bound the friction as in dynbfr) 182 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 183 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 184 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 185 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 186 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 187 END DO 188 END DO 189 189 190 190 ! Set surface roughness length 191 191 SELECT CASE ( nn_z0_met ) 192 192 ! 193 CASE ( 0 ) ! Constant roughness 193 CASE ( 0 ) ! Constant roughness 194 194 zhsro(:,:) = rn_hsro 195 195 CASE ( 1 ) ! Standard Charnock formula … … 227 227 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 228 228 DO jk = 2, jpkm1 229 DO jj = 2, jpjm1 229 DO jj = 2, jpjm1 230 230 DO ji = fs_2, fs_jpim1 ! vector opt. 231 231 zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1) … … 276 276 IF( ln_sigpsi ) THEN 277 277 zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1. 278 zwall_psi(ji,jj,jk) = rsc_psi / & 278 zwall_psi(ji,jj,jk) = rsc_psi / & 279 279 & ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 280 280 ELSE … … 295 295 ! diagonal 296 296 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & 297 & + rdt * zdiss * tmask(ji,jj,jk) 297 & + rdt * zdiss * tmask(ji,jj,jk) 298 298 ! 299 299 ! right hand side in en … … 317 317 ! First level 318 318 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 319 en(:,:,1) = MAX(en(:,:,1), rn_emin) 319 en(:,:,1) = MAX(en(:,:,1), rn_emin) 320 320 z_elem_a(:,:,1) = en(:,:,1) 321 321 z_elem_c(:,:,1) = 0._wp 322 322 z_elem_b(:,:,1) = 1._wp 323 ! 323 ! 324 324 ! One level below 325 325 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 326 326 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 327 327 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 328 z_elem_a(:,:,2) = 0._wp 328 z_elem_a(:,:,2) = 0._wp 329 329 z_elem_c(:,:,2) = 0._wp 330 330 z_elem_b(:,:,2) = 1._wp … … 335 335 ! Dirichlet conditions at k=1 336 336 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 337 en(:,:,1) = MAX(en(:,:,1), rn_emin) 337 en(:,:,1) = MAX(en(:,:,1), rn_emin) 338 338 z_elem_a(:,:,1) = en(:,:,1) 339 339 z_elem_c(:,:,1) = 0._wp … … 358 358 SELECT CASE ( nn_bc_bot ) 359 359 ! 360 CASE ( 0 ) ! Dirichlet 360 CASE ( 0 ) ! Dirichlet 361 361 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 362 362 ! ! Balance between the production and the dissipation terms … … 378 378 z_elem_c(ji,jj,ibotm1) = 0._wp 379 379 z_elem_b(ji,jj,ibotm1) = 1._wp 380 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 380 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 381 381 END DO 382 382 END DO 383 383 ! 384 384 CASE ( 1 ) ! Neumman boundary condition 385 ! 385 ! 386 386 !CDIR NOVERRCHK 387 387 DO jj = 2, jpjm1 … … 429 429 END DO 430 430 END DO 431 ! ! set the minimum value of tke 431 ! ! set the minimum value of tke 432 432 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 433 433 … … 471 471 DO jj = 2, jpjm1 472 472 DO ji = fs_2, fs_jpim1 ! vector opt. 473 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 473 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 474 474 END DO 475 475 END DO … … 490 490 ! 491 491 ! psi / k 492 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 492 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 493 493 ! 494 494 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable) … … 510 510 zesh2 = dir * ( prod + buoy ) + (1._wp - dir ) * prod ! production term 511 511 zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 512 ! 512 ! 513 513 ! building the matrix 514 514 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) … … 552 552 z_elem_c(:,:,2) = 0._wp 553 553 z_elem_b(:,:,2) = 1._wp 554 ! 554 ! 555 555 ! 556 556 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz … … 576 576 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 577 577 578 ! 578 ! 579 579 ! 580 580 END SELECT … … 586 586 ! 587 587 ! 588 CASE ( 0 ) ! Dirichlet 588 CASE ( 0 ) ! Dirichlet 589 589 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 590 590 ! ! Balance between the production and the dissipation terms … … 611 611 ! 612 612 CASE ( 1 ) ! Neumman boundary condition 613 ! 613 ! 614 614 !CDIR NOVERRCHK 615 615 DO jj = 2, jpjm1 … … 693 693 DO jj = 2, jpjm1 694 694 DO ji = fs_2, fs_jpim1 ! vector opt. 695 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 695 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 696 696 END DO 697 697 END DO … … 720 720 eps(ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 721 721 mxln(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 722 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 722 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 723 723 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 724 724 IF (ln_length_lim) mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) ) 725 725 END DO 726 726 END DO 727 END DO 727 END DO 728 728 729 729 ! … … 807 807 END DO 808 808 END DO 809 #if defined key_traldf_c2d || key_traldf_c3d 809 810 IF( ln_stopack) THEN 810 811 IF( nn_spp_avt > 0 ) CALL spp_gen(kt, avt(:,:,jk), nn_spp_avt, rn_avt_sd,jk) 811 812 IF( nn_spp_avm > 0 ) CALL spp_gen(kt, avm(:,:,jk), nn_spp_avm, rn_avm_sd,jk) 812 ENDIF 813 ENDIF 814 #else 815 IF ( ln_stopack ) & 816 & CALL ctl_stop( 'zdf_gls: parameter perturbation will only work with '// & 817 'key_traldf_c2d or key_traldf_c3d') 818 #endif 813 819 END DO 814 820 ! … … 851 857 !!---------------------------------------------------------------------- 852 858 !! *** ROUTINE zdf_gls_init *** 853 !! 854 !! ** Purpose : Initialization of the vertical eddy diffivity and 859 !! 860 !! ** Purpose : Initialization of the vertical eddy diffivity and 855 861 !! viscosity when using a gls turbulent closure scheme 856 862 !! … … 1071 1077 ! 1072 1078 END SELECT 1073 1079 1074 1080 ! !* Set Schmidt number for psi diffusion in the wave breaking case 1075 1081 ! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 1076 1082 ! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001 1077 1083 IF( ln_sigpsi ) THEN 1078 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1084 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1079 1085 ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 1080 1086 ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work … … 1084 1090 rsc_psi0 = rsc_psi 1085 1091 ENDIF 1086 1092 1087 1093 ! !* Shear free turbulence parameters 1088 1094 ! … … 1130 1136 rc04 = rc03 * rc0 1131 1137 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking 1132 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking 1138 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking 1133 1139 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1134 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1140 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1135 1141 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1136 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1142 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1137 1143 rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 1138 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1144 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1139 1145 1140 1146 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke … … 1151 1157 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 1152 1158 END DO 1153 ! 1159 ! 1154 1160 CALL gls_rst( nit000, 'READ' ) !* read or initialize all required files 1155 1161 ! … … 1162 1168 !!--------------------------------------------------------------------- 1163 1169 !! *** ROUTINE ts_rst *** 1164 !! 1170 !! 1165 1171 !! ** Purpose : Read or write TKE file (en) in restart file 1166 1172 !! 1167 1173 !! ** Method : use of IOM library 1168 !! if the restart does not contain TKE, en is either 1174 !! if the restart does not contain TKE, en is either 1169 1175 !! set to rn_emin or recomputed (nn_igls/=0) 1170 1176 !!---------------------------------------------------------------------- … … 1178 1184 !!---------------------------------------------------------------------- 1179 1185 ! 1180 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 1186 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 1181 1187 ! ! --------------- 1182 1188 IF( ln_rstart ) THEN !* Read the restart file … … 1197 1203 CALL iom_get( numror, jpdom_autoglo, 'mxln' , mxln ) 1198 1204 IF(nn_timing == 2) CALL timing_stop('iom_rstget') 1199 ELSE 1205 ELSE 1200 1206 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 1201 1207 en (:,:,:) = rn_emin 1202 mxln(:,:,:) = 0.05 1208 mxln(:,:,:) = 0.05 1203 1209 avt_k (:,:,:) = avt (:,:,:) 1204 1210 avm_k (:,:,:) = avm (:,:,:) … … 1210 1216 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' 1211 1217 en (:,:,:) = rn_emin 1212 mxln(:,:,:) = 0.05 1218 mxln(:,:,:) = 0.05 1213 1219 ENDIF 1214 1220 ! … … 1217 1223 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1218 1224 IF(nn_timing == 2) CALL timing_start('iom_rstput') 1219 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1225 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1220 1226 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 1221 1227 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 1222 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1228 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1223 1229 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 1224 1230 CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln ) … … 1239 1245 WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?' 1240 1246 END SUBROUTINE zdf_gls_init 1241 1247 1242 1248 SUBROUTINE zdf_gls( kt ) ! Empty routine 1243 1249 IMPLICIT NONE 1244 INTEGER, INTENT(in) :: kt ! ocean time step 1250 INTEGER, INTENT(in) :: kt ! ocean time step 1245 1251 WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt 1246 1252 END SUBROUTINE zdf_gls 1247 1253 1248 1254 SUBROUTINE gls_rst( kt, cdrw ) ! Empty routine 1249 1255 IMPLICIT NONE … … 1256 1262 !!====================================================================== 1257 1263 END MODULE zdfgls 1258 -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r11442 r13355 2 2 !!====================================================================== 3 3 !! *** MODULE zdftke *** 4 !! Ocean physics: vertical mixing coefficient computed from the tke 4 !! Ocean physics: vertical mixing coefficient computed from the tke 5 5 !! turbulent closure parameterization 6 6 !!===================================================================== … … 22 22 !! - ! 2008-05 (J.-M. Molines, G. Madec) 2D form of avtb 23 23 !! - ! 2008-06 (G. Madec) style + DOCTOR name for namelist parameters 24 !! - ! 2008-12 (G. Reffray) stable discretization of the production term 25 !! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl 24 !! - ! 2008-12 (G. Reffray) stable discretization of the production term 25 !! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl 26 26 !! ! + cleaning of the parameters + bugs correction 27 27 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase … … 52 52 USE wrk_nemo ! work arrays 53 53 USE timing ! Timing 54 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 54 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 55 55 #if defined key_agrif 56 56 USE agrif_opa_interp … … 74 74 INTEGER :: nn_pdl ! Prandtl number or not (ratio avt/avm) (=0/1) 75 75 REAL(wp) :: rn_ediff ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 76 REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation 76 REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation 77 77 REAL(wp) :: rn_ebb ! coefficient of the surface input of tke 78 78 REAL(wp) :: rn_emin ! minimum value of tke [m2/s2] … … 102 102 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_niw !: TKE budget- near-inertial waves term 103 103 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: efr ! surface boundary condition for nn_etau = 4 104 104 105 105 !! * Substitutions 106 106 # include "domzgr_substitute.h90" … … 118 118 !!---------------------------------------------------------------------- 119 119 ALLOCATE( & 120 & efr (jpi,jpj) , e_niw(jpi,jpj,jpk) , & 120 & efr (jpi,jpj) , e_niw(jpi,jpj,jpk) , & 121 121 #if defined key_c1d 122 122 & e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , & … … 147 147 !! surface: en = max( rn_emin0, rn_ebb * taum ) 148 148 !! bottom : en = rn_emin 149 !! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 150 !! 151 !! The now Turbulent kinetic energy is computed using the following 149 !! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 150 !! 151 !! The now Turbulent kinetic energy is computed using the following 152 152 !! time stepping: implicit for vertical diffusion term, linearized semi 153 !! implicit for kolmogoroff dissipation term, and explicit forward for 154 !! both buoyancy and shear production terms. Therefore a tridiagonal 153 !! implicit for kolmogoroff dissipation term, and explicit forward for 154 !! both buoyancy and shear production terms. Therefore a tridiagonal 155 155 !! linear system is solved. Note that buoyancy and shear terms are 156 156 !! discretized in a energy conserving form (Bruchard 2002). … … 160 160 !! 161 161 !! The now vertical eddy vicosity and diffusivity coefficients are 162 !! given by: 162 !! given by: 163 163 !! avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 164 !! avt = max( avmb, pdl * avm ) 164 !! avt = max( avmb, pdl * avm ) 165 165 !! eav = max( avmb, avm ) 166 166 !! where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and 167 !! given by an empirical funtion of the localRichardson number if nn_pdl=1 167 !! given by an empirical funtion of the localRichardson number if nn_pdl=1 168 168 !! 169 169 !! ** Action : compute en (now turbulent kinetic energy) … … 180 180 ! 181 181 IF( kt /= nit000 ) THEN ! restore before value to compute tke 182 avt (:,:,:) = avt_k (:,:,:) 183 avm (:,:,:) = avm_k (:,:,:) 184 avmu(:,:,:) = avmu_k(:,:,:) 185 avmv(:,:,:) = avmv_k(:,:,:) 186 ENDIF 187 ! 182 avt (:,:,:) = avt_k (:,:,:) 183 avm (:,:,:) = avm_k (:,:,:) 184 avmu(:,:,:) = avmu_k(:,:,:) 185 avmv(:,:,:) = avmv_k(:,:,:) 186 ENDIF 187 ! 188 #if defined key_traldf_c2d || key_traldf_c3d 188 189 IF( ln_stopack) THEN 189 190 IF( nn_spp_tkelc > 0 ) THEN … … 208 209 ENDIF 209 210 ENDIF 211 #else 212 IF ( ln_stopack ) & 213 & CALL ctl_stop( 'zdf_tke: parameter perturbation will only work with '// & 214 'key_traldf_c2d or key_traldf_c3d') 215 #endif 210 216 ! 211 217 CALL tke_tke ! now tke (en) … … 213 219 CALL tke_avn ! now avt, avm, avmu, avmv 214 220 ! 215 avt_k (:,:,:) = avt (:,:,:) 216 avm_k (:,:,:) = avm (:,:,:) 217 avmu_k(:,:,:) = avmu(:,:,:) 218 avmv_k(:,:,:) = avmv(:,:,:) 221 avt_k (:,:,:) = avt (:,:,:) 222 avm_k (:,:,:) = avm (:,:,:) 223 avmu_k(:,:,:) = avmu(:,:,:) 224 avmv_k(:,:,:) = avmv(:,:,:) 219 225 ! 220 226 #if defined key_agrif 221 ! Update child grid f => parent grid 227 ! Update child grid f => parent grid 222 228 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 223 #endif 229 #endif 224 230 IF ( kt == nitend ) THEN 225 231 DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 ) … … 271 277 CALL wrk_alloc( jpi,jpj, zfact2 ) 272 278 CALL wrk_alloc( jpi,jpj, zfact3 ) 273 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 279 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 274 280 ! 275 281 zbbrau = rn_ebb0 / rau0 ! Local constant initialisation 276 zfact1 = -.5_wp * rdt 282 zfact1 = -.5_wp * rdt 277 283 zfact2 = 1.5_wp * rdt * rn_ediss0 278 284 zfact3 = 0.5_wp * rn_ediss0 … … 294 300 END DO 295 301 END DO 296 302 297 303 !!bfr - start commented area 298 304 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 333 339 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 334 340 DO jk = jpkm1, 2, -1 335 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 341 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 336 342 DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1) 337 343 zus = zcof * taum(ji,jj) … … 341 347 END DO 342 348 ! ! finite LC depth 343 DO jj = 1, jpj 349 DO jj = 1, jpj 344 350 DO ji = 1, jpi 345 351 zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) … … 378 384 DO ji = 1, jpi 379 385 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 380 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 386 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 381 387 & / ( fse3uw_n(ji,jj,jk) & 382 388 & * fse3uw_b(ji,jj,jk) ) … … 407 413 ! ! shear prod. at w-point weightened by mask 408 414 zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 409 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 415 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 410 416 ! 411 417 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) … … 464 470 END DO 465 471 466 ! ! Save TKE prior to nn_etau addition 467 e_niw(:,:,:) = en(:,:,:) 468 ! 472 ! ! Save TKE prior to nn_etau addition 473 e_niw(:,:,:) = en(:,:,:) 474 ! 469 475 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 470 476 ! ! TKE due to surface and internal wave breaking … … 508 514 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 509 515 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 510 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 511 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 516 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 517 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 512 518 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 513 519 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & … … 517 523 END DO 518 524 ELSEIF( nn_etau == 4 ) THEN !* column integral independant of htau (rn_efr must be scaled up) 519 IF( nn_htau == 2 ) THEN ! efr dependant on time-varying htau 525 IF( nn_htau == 2 ) THEN ! efr dependant on time-varying htau 520 526 DO jj = 2, jpjm1 521 527 DO ji = fs_2, fs_jpim1 ! vector opt. … … 535 541 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 536 542 ! 537 DO jk = 2, jpkm1 ! TKE budget: near-inertial waves term 538 DO jj = 2, jpjm1 539 DO ji = fs_2, fs_jpim1 ! vector opt. 540 e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk) 541 END DO 542 END DO 543 END DO 544 ! 545 CALL lbc_lnk( e_niw, 'W', 1. ) 543 DO jk = 2, jpkm1 ! TKE budget: near-inertial waves term 544 DO jj = 2, jpjm1 545 DO ji = fs_2, fs_jpim1 ! vector opt. 546 e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk) 547 END DO 548 END DO 549 END DO 550 ! 551 CALL lbc_lnk( e_niw, 'W', 1. ) 546 552 ! 547 553 CALL wrk_dealloc( jpi,jpj, imlc ) ! integer 548 CALL wrk_dealloc( jpi,jpj, zhlc ) 549 CALL wrk_dealloc( jpi,jpj, zbbrau ) 550 CALL wrk_dealloc( jpi,jpj, zfact2 ) 551 CALL wrk_dealloc( jpi,jpj, zfact3 ) 552 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 554 CALL wrk_dealloc( jpi,jpj, zhlc ) 555 CALL wrk_dealloc( jpi,jpj, zbbrau ) 556 CALL wrk_dealloc( jpi,jpj, zfact2 ) 557 CALL wrk_dealloc( jpi,jpj, zfact3 ) 558 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 553 559 ! 554 560 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') … … 563 569 !! ** Purpose : Compute the vertical eddy viscosity and diffusivity 564 570 !! 565 !! ** Method : At this stage, en, the now TKE, is known (computed in 566 !! the tke_tke routine). First, the now mixing lenth is 571 !! ** Method : At this stage, en, the now TKE, is known (computed in 572 !! the tke_tke routine). First, the now mixing lenth is 567 573 !! computed from en and the strafification (N^2), then the mixings 568 574 !! coefficients are computed. 569 575 !! - Mixing length : a first evaluation of the mixing lengh 570 576 !! scales is: 571 !! mxl = sqrt(2*en) / N 577 !! mxl = sqrt(2*en) / N 572 578 !! where N is the brunt-vaisala frequency, with a minimum value set 573 579 !! to rmxl_min (rn_mxl0) in the interior (surface) ocean. 574 !! The mixing and dissipative length scale are bound as follow : 580 !! The mixing and dissipative length scale are bound as follow : 575 581 !! nn_mxl=0 : mxl bounded by the distance to surface and bottom. 576 582 !! zmxld = zmxlm = mxl 577 583 !! nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl 578 !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 584 !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 579 585 !! less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl 580 586 !! nn_mxl=3 : mxl is bounded from the surface to the bottom usings 581 !! |d/dz(xml)|<1 to obtain lup, and from the bottom to 582 !! the surface to obtain ldown. the resulting length 587 !! |d/dz(xml)|<1 to obtain lup, and from the bottom to 588 !! the surface to obtain ldown. the resulting length 583 589 !! scales are: 584 !! zmxld = sqrt( lup * ldown ) 590 !! zmxld = sqrt( lup * ldown ) 585 591 !! zmxlm = min ( lup , ldown ) 586 592 !! - Vertical eddy viscosity and diffusivity: 587 593 !! avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 588 !! avt = max( avmb, pdlr * avm ) 594 !! avt = max( avmb, pdlr * avm ) 589 595 !! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. 590 596 !! … … 601 607 IF( nn_timing == 1 ) CALL timing_start('tke_avn') 602 608 603 CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 609 CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 604 610 605 611 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 610 616 ! 611 617 ! initialisation of interior minimum value (avoid a 2d loop with mikt) 612 zmxlm(:,:,:) = rmxl_min 618 zmxlm(:,:,:) = rmxl_min 613 619 zmxld(:,:,:) = rmxl_min 614 620 ! … … 620 626 END DO 621 627 END DO 622 ELSE 628 ELSE 623 629 zmxlm(:,:,1) = rn_mxl0 624 630 ENDIF … … 638 644 ! !* Physical limits for the mixing length 639 645 ! 640 zmxld(:,:,1 ) = zmxlm(:,:,1) ! surface set to the minimum value 646 zmxld(:,:,1 ) = zmxlm(:,:,1) ! surface set to the minimum value 641 647 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 642 648 ! … … 738 744 END DO 739 745 END DO 746 #if defined key_traldf_c2d || key_traldf_c3d 740 747 IF( ln_stopack) THEN 741 748 IF(nn_spp_avt > 0 ) CALL spp_gen( 1, avt(:,:,jk), nn_spp_avt, rn_avt_sd, jk_spp_avt, jk) 742 749 IF(nn_spp_avm > 0 ) CALL spp_gen( 1, avm(:,:,jk), nn_spp_avm, rn_avm_sd, jk_spp_avm, jk) 743 750 ENDIF 751 #else 752 IF ( ln_stopack ) & 753 & CALL ctl_stop( 'tke_avn: parameter perturbation will only work with '// & 754 'key_traldf_c2d or key_traldf_c3d') 755 #endif 744 756 END DO 745 757 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) … … 787 799 ENDIF 788 800 ! 789 CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 801 CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 790 802 ! 791 803 IF( nn_timing == 1 ) CALL timing_stop('tke_avn') … … 797 809 !!---------------------------------------------------------------------- 798 810 !! *** ROUTINE zdf_tke_init *** 799 !! 800 !! ** Purpose : Initialization of the vertical eddy diffivity and 811 !! 812 !! ** Purpose : Initialization of the vertical eddy diffivity and 801 813 !! viscosity when using a tke turbulent closure scheme 802 814 !! … … 814 826 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 815 827 & rn_mxl0 , nn_pdl , ln_lc , rn_lc , & 816 & nn_etau , nn_htau , rn_efr , rn_c 828 & nn_etau , nn_htau , rn_efr , rn_c 817 829 !!---------------------------------------------------------------------- 818 830 … … 884 896 ALLOCATE( rn_ebb0 (jpi,jpj) ) ; rn_ebb0 = rn_ebb 885 897 ALLOCATE( rn_efr0 (jpi,jpj) ) ; rn_efr0 = rn_efr 886 898 887 899 IF( nn_etau == 2 ) THEN 888 900 ierr = zdf_mxl_alloc() … … 896 908 897 909 ! !* depth of penetration of surface tke 898 IF( nn_etau /= 0 ) THEN 910 IF( nn_etau /= 0 ) THEN 899 911 htau(:,:) = 0._wp 900 912 SELECT CASE( nn_htau ) ! Choice of the depth of penetration … … 902 914 htau(:,:) = 10._wp 903 915 CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees 904 htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 916 htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 905 917 CASE( 2 ) ! fraction of depth-integrated TKE within mixed-layer 906 918 rhtau = -1._wp / LOG( 1._wp - rn_c ) … … 945 957 END DO 946 958 dissl(:,:,:) = 1.e-12_wp 947 ! 959 ! 948 960 CALL tke_rst( nit000, 'READ' ) !* read or initialize all required files 949 961 ! … … 954 966 !!--------------------------------------------------------------------- 955 967 !! *** ROUTINE tke_rst *** 956 !! 968 !! 957 969 !! ** Purpose : Read or write TKE file (en) in restart file 958 970 !! 959 971 !! ** Method : use of IOM library 960 !! if the restart does not contain TKE, en is either 961 !! set to rn_emin or recomputed 972 !! if the restart does not contain TKE, en is either 973 !! set to rn_emin or recomputed 962 974 !!---------------------------------------------------------------------- 963 975 INTEGER , INTENT(in) :: kt ! ocean time-step … … 968 980 !!---------------------------------------------------------------------- 969 981 ! 970 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 982 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 971 983 ! ! --------------- 972 984 IF( ln_rstart ) THEN !* Read the restart file … … 1006 1018 avmu_k(:,:,:) = avmu(:,:,:) 1007 1019 avmv_k(:,:,:) = avmv(:,:,:) 1008 1020 1009 1021 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 1010 1022 ! ! ------------------- -
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/step.F90
r11442 r13355 110 110 IF( lk_tide ) CALL sbc_tide( kstp ) 111 111 IF( lk_bdy ) THEN 112 IF( ln_apr_dyn) CALL sbc_apr( kstp ) ! bdy_dta needs ssh_ib 112 IF( ln_apr_dyn) CALL sbc_apr( kstp ) ! bdy_dta needs ssh_ib 113 113 CALL bdy_dta ( kstp, time_offset=+1 ) ! update dynamic & tracer data at open boundaries 114 114 ENDIF … … 118 118 CALL lbc_lnk( tsb(:,:,:,tind), 'T', 1. ) 119 119 END DO 120 120 121 121 IF( ln_stopack ) CALL stopack_pert( kstp ) 122 122 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) … … 153 153 ENDIF 154 154 IF( ln_rnf_mouth ) THEN ! increase diffusivity at rivers mouths 155 #if defined key_traldf_c2d || key_traldf_c3d 155 156 IF ( ln_stopack .AND. ( nn_spp_arnf > 0 ) ) THEN 156 157 rn_avt_rnf0 = rn_avt_rnf 157 158 CALL spp_gen( kstp, rn_avt_rnf0,nn_spp_arnf,rn_arnf_sd,jk_spp_arnf ) 158 159 ENDIF 160 #else 161 IF ( ln_stopack .AND. ( nn_spp_arnf > 0 ) ) & 162 & CALL ctl_stop( 'stp: parameter perturbation will only work with '// & 163 'key_traldf_c2d or key_traldf_c3d') 164 #endif 159 165 DO jk = 2, nkrnf ; avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf0(:,:) * rnfmsk(:,:) * tmask(:,:,jk) ; END DO 160 166 ENDIF … … 204 210 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 205 211 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_cur) 206 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 207 CALL wzv ( kstp ) ! now cross-level velocity 208 209 IF( lk_dynspg_ts ) THEN 212 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 213 CALL wzv ( kstp ) ! now cross-level velocity 214 215 IF( lk_dynspg_ts ) THEN 210 216 ! In case the time splitting case, update almost all momentum trends here: 211 217 ! Note that the computation of vertical velocity above, hence "after" sea level … … 244 250 CALL div_cur( kstp ) ! Horizontal divergence & Relative vorticity (2nd call in time-split case) 245 251 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 246 CALL wzv ( kstp ) ! now cross-level velocity 252 CALL wzv ( kstp ) ! now cross-level velocity 247 253 ENDIF 248 254 … … 320 326 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 321 327 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 322 IF( ln_zps .AND. ln_isfcav) & 328 IF( ln_zps .AND. ln_isfcav) & 323 329 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps for top cell (ISF) 324 330 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & … … 334 340 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 335 341 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 336 IF( ln_zps .AND. ln_isfcav) & 342 IF( ln_zps .AND. ln_isfcav) & 337 343 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps for top cell (ISF) 338 344 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & … … 349 355 ua(:,:,:) = ua_sv(:,:,:) 350 356 va(:,:,:) = va_sv(:,:,:) 351 ! Revert now divergence and rotational to previously computed ones 357 ! Revert now divergence and rotational to previously computed ones 352 358 !(needed because of the time swap in div_cur, at the beginning of each time step) 353 359 hdivn(:,:,:) = hdivb(:,:,:) 354 rotn(:,:,:) = rotb(:,:,:) 360 rotn(:,:,:) = rotb(:,:,:) 355 361 356 362 CALL dyn_bfr( kstp ) ! bottom friction … … 398 404 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 399 405 ! AGRIF 400 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 401 CALL Agrif_Integrate_ChildGrids( stp ) 406 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 407 CALL Agrif_Integrate_ChildGrids( stp ) 402 408 403 409 IF ( Agrif_NbStepint().EQ.0 ) THEN … … 431 437 ! 432 438 #if defined key_iomput 433 IF( kstp == nitend .OR. indic < 0 ) THEN 439 IF( kstp == nitend .OR. indic < 0 ) THEN 434 440 CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF 435 IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 441 IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 436 442 ENDIF 437 443 #endif 438 444 ! 439 445 IF( nn_timing == 1 .AND. kstp == nit000 ) CALL timing_reset 440 ! 446 ! 441 447 ! 442 448 END SUBROUTINE stp
Note: See TracChangeset
for help on using the changeset viewer.