- Timestamp:
- 2017-12-12T16:38:41+01:00 (7 years ago)
- Location:
- branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 34 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r8552 r8992 20 20 !! dyn_asm_inc : Apply the dynamic (u and v) increments 21 21 !! ssh_asm_inc : Apply the SSH increment 22 !! ssh_asm_div : Apply divergence associated with SSH increment 22 23 !! seaice_asm_inc : Apply the seaice increment 23 24 !!---------------------------------------------------------------------- … … 48 49 PUBLIC dyn_asm_inc !: Apply the dynamic (u and v) increments 49 50 PUBLIC ssh_asm_inc !: Apply the SSH increment 51 PUBLIC ssh_asm_div !: Apply the SSH divergence 50 52 PUBLIC seaice_asm_inc !: Apply the seaice increment 51 53 … … 798 800 END SUBROUTINE ssh_asm_inc 799 801 802 SUBROUTINE ssh_asm_div( kt, phdivn ) 803 !!---------------------------------------------------------------------- 804 !! *** ROUTINE ssh_asm_div *** 805 !! 806 !! ** Purpose : ssh increment with z* is incorporated via a correction of the local divergence 807 !! across all the water column 808 !! 809 !! ** Method : 810 !! CAUTION : sshiau is positive (inflow) decreasing the 811 !! divergence and expressed in m/s 812 !! 813 !! ** Action : phdivn decreased by the ssh increment 814 !!---------------------------------------------------------------------- 815 INTEGER, INTENT(IN) :: kt ! ocean time-step index 816 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence 817 !! 818 INTEGER :: jk ! dummy loop index 819 REAL(wp), DIMENSION(:,:) , POINTER :: ztim ! local array 820 !!---------------------------------------------------------------------- 821 ! 822 #if defined key_asminc 823 CALL ssh_asm_inc( kt ) !== (calculate increments) 824 ! 825 IF( ln_linssh ) THEN 826 phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t_n(:,:,1) * tmask(:,:,1) 827 ELSE 828 CALL wrk_alloc( jpi,jpj, ztim) 829 ztim(:,:) = ssh_iau(:,:) / ( ht_n(:,:) + 1.0 - ssmask(:,:) ) 830 DO jk = 1, jpkm1 831 phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 832 END DO 833 ! 834 CALL wrk_dealloc( jpi,jpj, ztim) 835 ENDIF 836 #endif 837 ! 838 END SUBROUTINE ssh_asm_div 800 839 801 840 SUBROUTINE seaice_asm_inc( kt, kindic ) -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90
r7646 r8992 21 21 USE phycst ! physical constants 22 22 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 23 USE wet_dry ! Use wet dry to get reference ssh level 23 24 USE in_out_manager ! 24 25 … … 169 170 ii = idx%nbi(jb,igrd) 170 171 ij = idx%nbj(jb,igrd) 171 spgu(ii, ij) = dta%ssh(jb) 172 IF( ll_wd ) THEN 173 spgu(ii, ij) = dta%ssh(jb) - ssh_ref 174 ELSE 175 spgu(ii, ij) = dta%ssh(jb) 176 ENDIF 172 177 END DO 173 178 -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r8465 r8992 154 154 CALL iom_put( "e3tdef" , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 155 155 156 CALL iom_put( "ssh" , sshn ) ! sea surface height 156 IF( ll_wd ) THEN 157 CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying) 158 ELSE 159 CALL iom_put( "ssh" , sshn ) ! sea surface height 160 ENDIF 161 157 162 IF( iom_use("wetdep") ) & ! wet depth 158 CALL iom_put( "wetdep" , ht_ wd(:,:) + sshn(:,:) )163 CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) ) 159 164 160 165 CALL iom_put( "toce", tsn(:,:,:,jp_tem) ) ! 3D temperature -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r7822 r8992 39 39 USE domc1d ! 1D configuration: column location 40 40 USE dyncor_c1d ! 1D configuration: Coriolis term (cor_c1d routine) 41 USE wet_dry ! wetting and drying41 USE wet_dry, ONLY : ll_wd 42 42 ! 43 43 USE in_out_manager ! I/O manager … … 674 674 ENDIF 675 675 ! 676 IF( l n_wd ) THEN ! wetting and drying domain676 IF( ll_wd ) THEN ! wetting and drying domain 677 677 CALL iom_rstput( 0, 0, inum, 'ht_0' , ht_0 , ktype = jp_r8 ) 678 CALL iom_rstput( 0, 0, inum, 'ht_wd' , ht_wd , ktype = jp_r8 )679 678 ENDIF 680 679 ! -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r7753 r8992 680 680 ! 681 681 INTEGER :: ji, jj, jk ! dummy loop indices 682 REAL(wp) :: zlnwd ! =1./0. when ln_wd = T/F682 REAL(wp) :: zlnwd ! =1./0. when ln_wd_il = T/F 683 683 !!---------------------------------------------------------------------- 684 684 ! 685 685 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 686 686 ! 687 IF(ln_wd ) THEN687 IF(ln_wd_il) THEN 688 688 zlnwd = 1.0_wp 689 689 ELSE … … 876 876 ELSE !* Initialize at "rest" 877 877 ! 878 IF( ln_wd .AND. ( cn_cfg == 'wad' ) ) THEN 879 ! Wetting and drying test case 880 CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb ) 881 tsn (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones 882 sshn (:,:) = sshb(:,:) 883 un (:,:,:) = ub (:,:,:) 884 vn (:,:,:) = vb (:,:,:) 885 ! uniform T-S fields and initial ssh slope 886 ! needs to be called here and in istate which is called later. 887 ! Adjust vertical metrics 878 879 IF( ll_wd ) THEN ! MJB ll_wd edits start here - these are essential 880 ! 881 IF( cn_cfg == 'wad' ) THEN 882 ! Wetting and drying test case 883 CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb ) 884 tsn (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones 885 sshn (:,:) = sshb(:,:) 886 un (:,:,:) = ub (:,:,:) 887 vn (:,:,:) = vb (:,:,:) 888 ELSE 889 ! if not test case 890 sshn(:,:) = -ssh_ref 891 sshb(:,:) = -ssh_ref 892 893 DO jj = 1, jpj 894 DO ji = 1, jpi 895 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 896 897 sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 898 sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 899 ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 900 ENDIF 901 ENDDO 902 ENDDO 903 ENDIF !If test case else 904 905 ! Adjust vertical metrics for all wad 888 906 DO jk = 1, jpk 889 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &907 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 890 908 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 891 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))909 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 892 910 END DO 893 911 e3t_b(:,:,:) = e3t_n(:,:,:) 894 ! 895 ELSEIF( ln_wd ) THEN 896 ! 897 DO jj = 1, jpj 898 DO ji = 1, jpi 899 IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN 900 ! potential bug 901 ! Warning this assumes 2 layers only over wetting locations. needs investigating 902 e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 903 e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 904 e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 905 sshb(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) !!gm I don't understand that ! 906 sshn(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 907 ssha(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 908 ENDIF 909 ENDDO 910 ENDDO 912 913 DO ji = 1, jpi 914 DO jj = 1, jpj 915 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 916 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 917 ENDIF 918 END DO 919 END DO 920 911 921 ! 912 922 ELSE … … 916 926 sshn(:,:) = 0.0_wp 917 927 ! 918 END IF 928 END IF ! end of ll_wd edits 919 929 920 930 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r7646 r8992 18 18 USE dom_oce ! ocean space and time domain 19 19 USE phycst , ONLY : rsmall 20 USE wet_dry, ONLY : l n_wd, ht_wd20 USE wet_dry, ONLY : ll_wd ! Wetting and drying 21 21 ! 22 22 USE in_out_manager ! I/O manager … … 198 198 ENDIF 199 199 ! 200 IF( ln_wd ) THEN ! wetting and drying domain 201 CALL iom_rstput( 0, 0, inum, 'ht_0' , ht_0 , ktype = jp_r8 ) 202 CALL iom_rstput( 0, 0, inum, 'ht_wd' , ht_wd , ktype = jp_r8 ) 203 ENDIF 200 IF( ll_wd ) CALL iom_rstput( 0, 0, inum, 'ht_0' , ht_0 , ktype = jp_r8 ) 201 204 202 ! ! ============================ 205 203 CALL iom_close( inum ) ! close the files -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r7753 r8992 30 30 USE usrdef_zgr ! user defined vertical coordinate system 31 31 USE depth_e3 ! depth <=> e3 32 USE wet_dry, ONLY: ln_wd, ht_wd32 USE wet_dry, ONLY: ll_wd, ssh_ref ! Wetting and drying 33 33 ! 34 34 USE in_out_manager ! I/O manager … … 258 258 k_bot(:,:) = INT( z2d(:,:) ) 259 259 ! 260 ! bathymetry with orography (wetting and drying only)261 IF( l n_wd ) CALL iom_get( inum, jpdom_data, 'ht_wd' , ht_wd , lrowattr=ln_use_jattr)260 ! reference depth for negative bathy (wetting and drying only) 261 IF( ll_wd ) CALL iom_get( inum, 'rn_wd_ref_depth' , ssh_ref ) 262 262 ! 263 263 CALL iom_close( inum ) -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90
r7753 r8992 25 25 USE iscplhsb ! ice sheet / ocean coupling 26 26 USE iscplini ! ice sheet / ocean coupling 27 #if defined key_asminc 28 USE asminc ! Assimilation increment 29 #endif 27 30 ! 28 31 USE in_out_manager ! I/O manager … … 92 95 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) !== runoffs ==! (update hdivn field) 93 96 ! 97 #if defined key_asminc 98 IF( ln_sshinc .AND. ln_asmiau ) CALL ssh_asm_div( kt, hdivn ) !== SSH assimilation ==! (update hdivn field) 99 ! 100 #endif 94 101 IF( ln_isf ) CALL sbc_isf_div( hdivn ) !== ice shelf ==! (update hdivn field) 95 102 ! -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r7761 r8992 438 438 ! 439 439 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 440 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )440 IF( ln_wd_il ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 441 441 ! 442 442 IF( kt == nit000 ) THEN … … 451 451 ENDIF 452 452 ! 453 IF( ln_wd ) THEN453 IF( ln_wd_il ) THEN 454 454 DO jj = 2, jpjm1 455 455 DO ji = 2, jpim1 456 ll_tmp1 = MIN( sshn(ji,jj) ,sshn(ji+1,jj) ) > &457 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &458 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) )&459 & 460 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. (&461 & MAX( sshn(ji,jj) ,sshn(ji+1,jj) ) > &462 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )456 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 457 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 458 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 459 & > rn_wdmin1 + rn_wdmin2 460 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 461 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 462 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 463 463 464 464 IF(ll_tmp1) THEN … … 466 466 ELSE IF(ll_tmp2) THEN 467 467 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 468 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_ wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &469 & / (sshn(ji+1,jj) - 468 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 469 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 470 470 ELSE 471 471 zcpx(ji,jj) = 0._wp 472 472 END IF 473 473 474 ll_tmp1 = MIN( sshn(ji,jj) ,sshn(ji,jj+1) ) > &475 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &476 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) )&477 & 478 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. (&479 & MAX( sshn(ji,jj) ,sshn(ji,jj+1) ) > &480 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )474 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 475 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 476 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 477 & > rn_wdmin1 + rn_wdmin2 478 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 479 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 480 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 481 481 482 482 IF(ll_tmp1) THEN … … 484 484 ELSE IF(ll_tmp2) THEN 485 485 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 486 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_ wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &487 & / (sshn(ji,jj+1) - 486 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 487 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 488 488 ELSE 489 489 zcpy(ji,jj) = 0._wp … … 509 509 510 510 511 IF( ln_wd ) THEN511 IF( ln_wd_il ) THEN 512 512 513 513 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) … … 540 540 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 541 541 542 IF( ln_wd ) THEN542 IF( ln_wd_il ) THEN 543 543 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 544 544 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) … … 555 555 ! 556 556 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 557 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )557 IF( ln_wd_il ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 558 558 ! 559 559 END SUBROUTINE hpg_sco … … 700 700 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 701 701 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 702 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )703 ! 704 ! 705 IF( ln_wd ) THEN702 IF( ln_wd_il ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 703 ! 704 ! 705 IF( ln_wd_il ) THEN 706 706 DO jj = 2, jpjm1 707 707 DO ji = 2, jpim1 708 ll_tmp1 = MIN( sshn(ji,jj) ,sshn(ji+1,jj) ) > &709 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &710 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) )&711 & 712 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. (&713 & MAX( sshn(ji,jj) ,sshn(ji+1,jj) ) > &714 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )708 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 709 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 710 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 711 & > rn_wdmin1 + rn_wdmin2 712 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 713 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 714 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 715 715 716 716 IF(ll_tmp1) THEN … … 718 718 ELSE IF(ll_tmp2) THEN 719 719 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 720 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_ wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &721 & / (sshn(ji+1,jj) - 720 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 721 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 722 722 ELSE 723 723 zcpx(ji,jj) = 0._wp 724 724 END IF 725 725 726 ll_tmp1 = MIN( sshn(ji,jj) ,sshn(ji,jj+1) ) > &727 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &728 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) )&729 & 730 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. (&731 & MAX( sshn(ji,jj) ,sshn(ji,jj+1) ) > &732 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )726 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 727 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 728 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 729 & > rn_wdmin1 + rn_wdmin2 730 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 731 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 732 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 733 733 734 734 IF(ll_tmp1) THEN … … 736 736 ELSE IF(ll_tmp2) THEN 737 737 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 738 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_ wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &739 & / (sshn(ji,jj+1) - 738 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 739 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 740 740 ELSE 741 741 zcpy(ji,jj) = 0._wp … … 915 915 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 916 916 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 917 IF( ln_wd ) THEN917 IF( ln_wd_il ) THEN 918 918 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 919 919 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) … … 938 938 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 939 939 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 940 IF( ln_wd ) THEN940 IF( ln_wd_il ) THEN 941 941 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 942 942 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) … … 952 952 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 953 953 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 954 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )954 IF( ln_wd_il ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 955 955 ! 956 956 END SUBROUTINE hpg_djc … … 989 989 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 990 990 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 991 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )991 IF( ln_wd_il ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 992 992 ! 993 993 IF( kt == nit000 ) THEN … … 1002 1002 IF( ln_linssh ) znad = 0._wp 1003 1003 1004 IF( ln_wd ) THEN1004 IF( ln_wd_il ) THEN 1005 1005 DO jj = 2, jpjm1 1006 1006 DO ji = 2, jpim1 1007 ll_tmp1 = MIN( sshn(ji,jj) ,sshn(ji+1,jj) ) > &1008 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &1009 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) )&1010 & 1011 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. (&1012 & MAX( sshn(ji,jj) ,sshn(ji+1,jj) ) > &1013 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )1007 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1008 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1009 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 1010 & > rn_wdmin1 + rn_wdmin2 1011 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 1012 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1013 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1014 1014 1015 1015 IF(ll_tmp1) THEN … … 1017 1017 ELSE IF(ll_tmp2) THEN 1018 1018 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 1019 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_ wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &1019 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 1020 1020 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 1021 1022 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1021 1023 ELSE 1022 1024 zcpx(ji,jj) = 0._wp 1023 1025 END IF 1024 1026 1025 ll_tmp1 = MIN( sshn(ji,jj) ,sshn(ji,jj+1) ) > &1026 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &1027 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) )&1028 & 1029 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. (&1030 & MAX( sshn(ji,jj) ,sshn(ji,jj+1) ) > &1031 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )1027 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1028 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1029 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1030 & > rn_wdmin1 + rn_wdmin2 1031 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 1032 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1033 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1032 1034 1033 1035 IF(ll_tmp1) THEN … … 1035 1037 ELSE IF(ll_tmp2) THEN 1036 1038 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 1037 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 1038 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1039 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 1040 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1041 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 1042 1039 1043 ELSE 1040 1044 zcpy(ji,jj) = 0._wp … … 1227 1231 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1228 1232 ENDIF 1229 IF( ln_wd ) THEN1230 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1231 zdpdx2 = zdpdx2 * zcpx(ji,jj) 1233 IF( ln_wd_il ) THEN 1234 zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj) 1235 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1232 1236 ENDIF 1233 1237 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) … … 1286 1290 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1287 1291 ENDIF 1288 IF( ln_wd ) THEN1289 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1290 zdpdy2 = zdpdy2 * zcpy(ji,jj) 1292 IF( ln_wd_il ) THEN 1293 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1294 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1291 1295 ENDIF 1292 1296 … … 1301 1305 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1302 1306 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1303 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )1307 IF( ln_wd_il ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 1304 1308 ! 1305 1309 END SUBROUTINE hpg_prj -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r7753 r8992 28 28 USE dom_oce ! ocean space and time domain 29 29 USE sbc_oce ! Surface boundary condition: ocean fields 30 USE sbcrnf ! river runoffs 30 31 USE phycst ! physical constants 31 32 USE dynadv ! dynamics: vector invariant versus flux form … … 220 221 IF ( .NOT. ln_isf ) THEN ! if no ice shelf melting 221 222 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) & 222 & - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 223 & ) * tmask(:,:,1) 224 IF( ln_rnf_depth ) THEN 225 DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too, not sure for ice shelf case yet 226 DO jj = 1, jpj 227 DO ji = 1, jpi 228 IF( mikt(ji,jj) <= jk .and. jk <= nk_rnf(ji,jj) ) THEN 229 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * (-rnf_b(ji,jj) + rnf(ji,jj))*(e3t_n(ji,jj,jk)/h_rnf(ji,jj) )*tmask(ji,jj,jk) 230 ENDIF 231 ENDDO 232 ENDDO 233 ENDDO 234 ELSE 235 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 236 ENDIF 223 237 ELSE ! if ice shelf melting 224 238 DO jj = 1, jpj -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r7831 r8992 1 1 MODULE dynspg_ts 2 3 !! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out ! 4 2 5 !!====================================================================== 3 6 !! *** MODULE dynspg_ts *** … … 137 140 INTEGER, INTENT(in) :: kt ! ocean time-step index 138 141 ! 139 LOGICAL :: ll_fw_start 140 LOGICAL :: ll_init 142 LOGICAL :: ll_fw_start ! if true, forward integration 143 LOGICAL :: ll_init ! if true, special startup of 2d equations 141 144 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 142 INTEGER :: ji, jj, jk, jn 143 INTEGER :: ikbu, ikbv, noffset 144 INTEGER :: iktu, iktv 145 INTEGER :: ji, jj, jk, jn ! dummy loop indices 146 INTEGER :: ikbu, ikbv, noffset ! local integers 147 INTEGER :: iktu, iktv ! local integers 145 148 REAL(wp) :: zmdi 146 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf 149 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 150 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 148 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 151 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 149 152 REAL(wp) :: zu_spg, zv_spg ! - - 150 REAL(wp) :: zhura, zhvra ! - - 151 REAL(wp) :: za0, za1, za2, za3 ! - - 153 REAL(wp) :: zhura, zhvra ! - - 154 REAL(wp) :: za0, za1, za2, za3 ! - - 155 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. 156 157 INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point 158 152 159 ! 153 160 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e … … 158 165 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 159 166 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 167 REAL(wp), POINTER, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 168 REAL(wp), POINTER, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask 160 169 !!---------------------------------------------------------------------- 161 170 ! … … 169 178 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 170 179 CALL wrk_alloc( jpi,jpj, zhf ) 171 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 180 IF( ln_wd_il ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 181 IF( ln_wd_dl ) CALL wrk_alloc( jpi, jpj, ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2) 182 183 172 184 ! 173 185 zmdi=1.e+20 ! missing data indicator for masking … … 178 190 z1_2 = 0.5_wp 179 191 zraur = 1._wp / rau0 180 ! ! reciprocal of baroclinic time step 192 zwdramp = r_rn_wdmin1 ! simplest ramp 193 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 194 ! ! reciprocal of baroclinic time step 181 195 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 182 196 ELSE ; z2dt_bf = 2.0_wp * rdt … … 184 198 z1_2dt_b = 1.0_wp / z2dt_bf 185 199 ! 186 ll_init = ln_bt_av 200 ll_init = ln_bt_av ! if no time averaging, then no specific restart 187 201 ll_fw_start = .FALSE. 188 ! 202 ! ! time offset in steps for bdy data update 189 203 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro 190 204 ELSE ; noffset = 0 191 205 ENDIF 192 206 ! 193 IF( kt == nit000 ) THEN 207 IF( kt == nit000 ) THEN !* initialisation 194 208 ! 195 209 IF(lwp) WRITE(numout,*) … … 399 413 ! ! ---------------------------------------------------- 400 414 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 401 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters415 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 402 416 DO jj = 2, jpjm1 403 417 DO ji = 2, jpim1 404 ll_tmp1 = MIN( sshn(ji,jj) ,sshn(ji+1,jj) ) > &405 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &406 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) )&407 & 408 ll_tmp2 = ( ABS( sshn(ji+1,jj) -sshn(ji ,jj)) > 1.E-12 ).AND.( &409 & MAX( sshn(ji,jj) ,sshn(ji+1,jj) ) > &410 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )418 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 419 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 420 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 421 & > rn_wdmin1 + rn_wdmin2 422 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 423 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 424 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 411 425 412 426 IF(ll_tmp1) THEN … … 414 428 ELSE IF(ll_tmp2) THEN 415 429 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 416 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 417 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 430 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 431 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 432 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 433 418 434 ELSE 419 435 zcpx(ji,jj) = 0._wp 420 436 END IF 421 437 422 ll_tmp1 = MIN( sshn(ji,jj) ,sshn(ji,jj+1) ) > &423 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &424 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) )&425 & 426 ll_tmp2 = ( ABS( sshn(ji,jj) -sshn(ji,jj+1)) > 1.E-12 ).AND.( &427 & MAX( sshn(ji,jj) ,sshn(ji,jj+1) ) > &428 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )438 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 439 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 440 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 441 & > rn_wdmin1 + rn_wdmin2 442 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 443 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 444 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 429 445 430 446 IF(ll_tmp1) THEN … … 432 448 ELSE IF(ll_tmp2) THEN 433 449 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 434 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 435 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 450 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 451 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 452 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 453 436 454 ELSE 437 455 zcpy(ji,jj) = 0._wp … … 443 461 DO ji = 2, jpim1 444 462 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 445 & * r1_e1u(ji,jj) * zcpx(ji,jj) 463 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 446 464 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 447 & * r1_e2v(ji,jj) * zcpy(ji,jj) 465 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 466 448 467 END DO 449 468 END DO … … 468 487 END DO 469 488 ! 470 ! 489 ! ! Add bottom stress contribution from baroclinic velocities: 471 490 IF (ln_bt_fw) THEN 472 491 DO jj = 2, jpjm1 … … 490 509 ! 491 510 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 492 IF( ln_wd ) THEN493 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 494 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 511 IF( ln_wd_il ) THEN 512 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) * wdrampu(ji,jj) 513 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) * wdrampv(ji,jj) 495 514 ELSE 496 515 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) … … 627 646 vn_adv(:,:) = 0._wp 628 647 ! ! ==================== ! 648 649 IF (ln_wd_dl) THEN 650 zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary) 651 zvwdmask(:,:) = 0._wp ! 652 zuwdav2(:,:) = 0._wp 653 zvwdav2(:,:) = 0._wp 654 END IF 655 656 629 657 DO jn = 1, icycle ! sub-time-step loop ! 630 658 ! ! ==================== ! … … 654 682 ! Extrapolate Sea Level at step jit+0.5: 655 683 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 656 ! 684 685 ! set wetting & drying mask at tracer points for this barotropic sub-step 686 IF ( ln_wd_dl ) THEN 687 688 IF ( ln_wd_dl_rmp ) THEN 689 DO jj = 1, jpj 690 DO ji = 1, jpi ! vector opt. 691 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 692 ! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 693 ztwdmask(ji,jj) = 1._wp 694 ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 695 ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1)) ) 696 ELSE 697 ztwdmask(ji,jj) = 0._wp 698 END IF 699 END DO 700 END DO 701 ELSE 702 DO jj = 1, jpj 703 DO ji = 1, jpi ! vector opt. 704 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 705 ztwdmask(ji,jj) = 1._wp 706 ELSE 707 ztwdmask(ji,jj) = 0._wp 708 END IF 709 END DO 710 END DO 711 END IF 712 713 END IF 714 715 657 716 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 658 717 DO ji = 2, fs_jpim1 ! Vector opt. … … 706 765 ENDIF 707 766 #endif 708 IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 709 ! 767 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 768 769 IF ( ln_wd_dl ) THEN 770 771 772 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 773 774 DO jj = 1, jpjm1 775 DO ji = 1, jpim1 776 IF ( zwx(ji,jj) > 0.0 ) THEN 777 zuwdmask(ji, jj) = ztwdmask(ji ,jj) 778 ELSE 779 zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 780 END IF 781 zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 782 un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 783 784 IF ( zwy(ji,jj) > 0.0 ) THEN 785 zvwdmask(ji, jj) = ztwdmask(ji, jj ) 786 ELSE 787 zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 788 END IF 789 zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 790 vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 791 END DO 792 END DO 793 794 795 END IF 796 710 797 ! Sum over sub-time-steps to compute advective velocities 711 798 za2 = wgtbtp2(jn) 712 799 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 713 800 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 714 ! 801 802 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True) 803 IF ( ln_wd_dl_bc ) THEN 804 zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 805 zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 806 END IF 807 715 808 ! Set next sea level: 716 809 DO jj = 2, jpjm1 … … 766 859 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 767 860 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 768 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters861 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 769 862 DO jj = 2, jpjm1 770 863 DO ji = 2, jpim1 771 864 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 772 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &773 & MAX( zsshp2_e(ji,jj) + ht_ wd(ji,jj), zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) ) &865 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 866 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 774 867 & > rn_wdmin1 + rn_wdmin2 775 868 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( & 776 869 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 777 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )870 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 778 871 779 872 IF(ll_tmp1) THEN … … 781 874 ELSE IF(ll_tmp2) THEN 782 875 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 783 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &876 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 784 877 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 785 878 ELSE … … 788 881 789 882 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 790 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &791 & MAX( zsshp2_e(ji,jj) + ht_ wd(ji,jj), zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) ) &883 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 884 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 792 885 & > rn_wdmin1 + rn_wdmin2 793 886 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( & 794 887 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 795 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )888 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 796 889 797 890 IF(ll_tmp1) THEN … … 799 892 ELSE IF(ll_tmp2) THEN 800 893 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 801 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &894 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 802 895 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 803 896 ELSE … … 809 902 ! 810 903 ! Compute associated depths at U and V points: 811 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN 904 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 812 905 ! 813 906 DO jj = 2, jpjm1 … … 886 979 ! 887 980 ! Add bottom stresses: 888 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 889 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 981 !jth do implicitly instead 982 IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs 983 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 984 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 985 ENDIF 890 986 ! 891 987 ! Add top stresses: … … 895 991 ! Surface pressure trend: 896 992 897 IF( ln_wd ) THEN993 IF( ln_wd_il ) THEN 898 994 DO jj = 2, jpjm1 899 995 DO ji = 2, jpim1 … … 919 1015 ! 920 1016 ! Set next velocities: 921 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1017 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 922 1018 DO jj = 2, jpjm1 923 1019 DO ji = fs_2, fs_jpim1 ! vector opt. … … 933 1029 & + zv_frc(ji,jj) ) & 934 1030 & ) * ssvmask(ji,jj) 1031 1032 !jth implicit bottom friction: 1033 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 1034 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * bfrua(ji,jj) * hur_e(ji,jj)) 1035 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * bfrva(ji,jj) * hvr_e(ji,jj)) 1036 ENDIF 1037 935 1038 END DO 936 1039 END DO 937 1040 ! 938 ELSE 1041 ELSE !* Flux form 939 1042 DO jj = 2, jpjm1 940 1043 DO ji = fs_2, fs_jpim1 ! vector opt. 941 1044 942 IF( ln_wd ) THEN 943 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 944 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 945 ELSE 946 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 947 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 948 END IF 1045 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 1046 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 1047 949 1048 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 950 1049 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) … … 964 1063 END DO 965 1064 ENDIF 966 ! 1065 1066 967 1067 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 968 IF( ln_wd ) THEN 969 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 970 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 971 ELSE 972 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 973 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 974 END IF 1068 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 1069 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 975 1070 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 976 1071 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) … … 1005 1100 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 1006 1101 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 1007 ELSE ! Sum transports 1008 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 1009 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 1010 ENDIF 1011 ! ! Sum sea level 1102 ELSE ! Sum transports 1103 IF ( .NOT.ln_wd_dl ) THEN 1104 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 1105 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 1106 ELSE 1107 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:) 1108 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:) 1109 END IF 1110 ENDIF 1111 ! ! Sum sea level 1012 1112 ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 1113 1013 1114 ! ! ==================== ! 1014 1115 END DO ! end loop ! … … 1033 1134 vb2_b(:,:) = zwy(:,:) 1034 1135 ENDIF 1136 1137 1035 1138 ! 1036 1139 ! Update barotropic trend: … … 1062 1165 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1063 1166 ENDIF 1064 ! 1167 1168 1169 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 1065 1170 DO jk = 1, jpkm1 1066 ! Correct velocities:1067 1171 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 1068 1172 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 1069 !1070 1173 END DO 1071 ! 1174 1175 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 1176 DO jk = 1, jpkm1 1177 un(:,:,jk) = ( un_adv(:,:) + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)) ) * umask(:,:,jk) 1178 vn(:,:,jk) = ( vn_adv(:,:) + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)) ) * vmask(:,:,jk) 1179 END DO 1180 END IF 1181 1182 1072 1183 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 1073 1184 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current … … 1097 1208 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 1098 1209 CALL wrk_dealloc( jpi,jpj, zhf ) 1099 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 1210 IF( ln_wd_il ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 1211 IF( ln_wd_dl ) CALL wrk_dealloc( jpi, jpj, ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 ) 1100 1212 ! 1101 1213 IF ( ln_diatmb ) THEN … … 1271 1383 rdtbt = rdt / REAL( nn_baro , wp ) 1272 1384 zcmax = zcmax * rdtbt 1273 1385 ! Print results 1274 1386 IF(lwp) WRITE(numout,*) 1275 1387 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r7753 r8992 27 27 USE agrif_opa_interp 28 28 #endif 29 #if defined key_asminc30 USE asminc ! Assimilation increment31 #endif32 29 ! 33 30 USE in_out_manager ! I/O manager … … 94 91 ! ! After Sea Surface Height ! 95 92 ! !------------------------------! 96 IF(ln_wd ) THEN93 IF(ln_wd_il) THEN 97 94 CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 98 95 ENDIF … … 121 118 ENDIF 122 119 ENDIF 123 124 #if defined key_asminc125 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN ! Include the IAU weighted SSH increment126 CALL ssh_asm_inc( kt )127 ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)128 ENDIF129 #endif130 120 ! !------------------------------! 131 121 ! ! outputs ! -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
r7646 r8992 1 1 MODULE wet_dry 2 3 !! includes updates to namelist namwad for diagnostic outputs of ROMS wetting and drying 4 2 5 !!============================================================================== 3 6 !! *** MODULE wet_dry *** 4 7 !! Wetting and drying includes initialisation routine and routines to 5 8 !! compute and apply flux limiters and preserve water depth positivity 6 !! only effects if wetting/drying is on (ln_wd == .true.)9 !! only effects if wetting/drying is on (ln_wd_il == .true. or ln_wd_dl==.true. ) 7 10 !!============================================================================== 8 11 !! History : 3.6 ! 2014-09 ((H.Liu) Original code … … 31 34 !! --------------------------------------------------------------------- 32 35 33 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdmask !: u- and v- limiter 34 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht_wd !: wetting and drying t-pt depths 35 ! (can include negative depths) 36 37 LOGICAL, PUBLIC :: ln_wd !: Wetting/drying activation switch (T:on,F:off) 36 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdmask !: u- and v- limiter (can include negative depths) 37 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdramp, wdrampu, wdrampv !: for hpg limiting 38 39 LOGICAL, PUBLIC :: ln_wd_il !: Wetting/drying il activation switch (T:on,F:off) 40 LOGICAL, PUBLIC :: ln_wd_dl !: Wetting/drying dl activation switch (T:on,F:off) 41 REAL(wp), PUBLIC :: rn_wdmin0 !: depth at which wetting/drying starts 38 42 REAL(wp), PUBLIC :: rn_wdmin1 !: minimum water depth on dried cells 39 REAL(wp), PUBLIC :: r n_wdmin2 !: tolerrance of minimum water depth on dried cells40 REAL(wp), PUBLIC :: rn_wd ld !: land elevation below which wetting/drying41 !:will be considered43 REAL(wp), PUBLIC :: r_rn_wdmin1 !: 1/minimum water depth on dried cells 44 REAL(wp), PUBLIC :: rn_wdmin2 !: tolerance of minimum water depth on dried cells 45 REAL(wp), PUBLIC :: rn_wdld !: land elevation below which wetting/drying will be considered 42 46 INTEGER , PUBLIC :: nn_wdit !: maximum number of iteration for W/D limiter 47 LOGICAL, PUBLIC :: ln_wd_dl_bc !: DL scheme: True implies 3D velocities are set to the barotropic values at points 48 !: where the flow is from wet points on less than half the barotropic sub-steps 49 LOGICAL, PUBLIC :: ln_wd_dl_rmp !: use a ramp for the dl flux limiter between 2 rn_wdmin1 and rn_wdmin1 (rather than a cut-off at rn_wdmin1) 50 REAL(wp), PUBLIC :: ssh_ref !: height of z=0 with respect to the geoid; 51 52 LOGICAL, PUBLIC :: ll_wd !: Wetting/drying activation switch if either ln_wd_il or ln_wd_dl 43 53 44 54 PUBLIC wad_init ! initialisation routine called by step.F90 … … 58 68 !! ** input : - namwad namelist 59 69 !!---------------------------------------------------------------------- 60 NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 70 NAMELIST/namwad/ ln_wd_il, ln_wd_dl, rn_wdmin0, rn_wdmin1, rn_wdmin2, rn_wdld, & 71 & nn_wdit, ln_wd_dl_bc, ln_wd_dl_rmp 72 61 73 INTEGER :: ios ! Local integer output status for namelist read 62 74 INTEGER :: ierr ! Local integer status array allocation … … 78 90 WRITE(numout,*) '~~~~~~~~' 79 91 WRITE(numout,*) ' Namelist namwad' 80 WRITE(numout,*) ' Logical activation ln_wd = ', ln_wd 92 WRITE(numout,*) ' Logical for Iter Lim wd option ln_wd_il = ', ln_wd_il 93 WRITE(numout,*) ' Logical for Dir. Lim wd option ln_wd_dl = ', ln_wd_dl 94 WRITE(numout,*) ' Depth at which wet/drying starts rn_wdmin0 = ', rn_wdmin0 81 95 WRITE(numout,*) ' Minimum wet depth on dried cells rn_wdmin1 = ', rn_wdmin1 82 WRITE(numout,*) ' Tolerance of min wet depth rn_wdmin2= ', rn_wdmin296 WRITE(numout,*) ' Tolerance of min wet depth rn_wdmin2 = ', rn_wdmin2 83 97 WRITE(numout,*) ' land elevation threshold rn_wdld = ', rn_wdld 84 98 WRITE(numout,*) ' Max iteration for W/D limiter nn_wdit = ', nn_wdit 99 WRITE(numout,*) ' T => baroclinic u,v=0 at dry pts: ln_wd_dl_bc = ', ln_wd_dl_bc 100 WRITE(numout,*) ' use a ramp for rwd limiter: ln_wd_dl_rwd_rmp = ', ln_wd_dl_rmp 101 85 102 ENDIF 86 ! 87 IF(ln_wd) THEN 88 ALLOCATE( wdmask(jpi,jpj), ht_wd(jpi,jpj), STAT=ierr ) 103 IF( .NOT. ln_read_cfg ) THEN 104 WRITE(numout,*) ' No configuration file so seting ssh_ref to zero ' 105 ssh_ref=0.0 106 ENDIF 107 108 r_rn_wdmin1=1/rn_wdmin1 109 ll_wd = .FALSE. 110 IF(ln_wd_il .OR. ln_wd_dl) THEN 111 ll_wd = .TRUE. 112 ALLOCATE( wdmask(jpi,jpj), STAT=ierr ) 113 ALLOCATE( wdramp(jpi,jpj), wdrampu(jpi,jpj), wdrampv(jpi,jpj), STAT=ierr ) 89 114 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 90 115 ENDIF … … 122 147 IF( nn_timing == 1 ) CALL timing_start('wad_lmt') 123 148 124 IF(ln_wd) THEN 125 126 CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 127 CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv) 128 ! 129 130 !IF(lwp) WRITE(numout,*) 131 !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting' 132 133 jflag = 0 134 zdepwd = 50._wp !maximum depth on which that W/D could possibly happen 135 136 137 zflxp(:,:) = 0._wp 138 zflxn(:,:) = 0._wp 139 zflxu(:,:) = 0._wp 140 zflxv(:,:) = 0._wp 141 142 zwdlmtu(:,:) = 1._wp 143 zwdlmtv(:,:) = 1._wp 144 145 ! Horizontal Flux in u and v direction 146 DO jk = 1, jpkm1 147 DO jj = 1, jpj 148 DO ji = 1, jpi 149 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 150 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 151 END DO 152 END DO 153 END DO 154 155 zflxu(:,:) = zflxu(:,:) * e2u(:,:) 156 zflxv(:,:) = zflxv(:,:) * e1v(:,:) 157 158 wdmask(:,:) = 1 159 DO jj = 2, jpj 160 DO ji = 2, jpi 161 162 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE ! we don't care about land cells 163 IF( ht_wd(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry 164 165 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & 166 & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji, jj-1), 0._wp) 167 zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj), 0._wp) + & 168 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 169 170 zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 171 IF(zdep2 .le. 0._wp) THEN !add more safty, but not necessary 172 sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 173 IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp 174 IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 175 IF(zflxv(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = 0._wp 176 IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 177 wdmask(ji,jj) = 0._wp 178 END IF 179 ENDDO 180 END DO 149 150 CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 151 CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv) 152 ! 153 154 jflag = 0 155 zdepwd = 50._wp !maximum depth on which that W/D could possibly happen 156 157 158 zflxp(:,:) = 0._wp 159 zflxn(:,:) = 0._wp 160 zflxu(:,:) = 0._wp 161 zflxv(:,:) = 0._wp 162 163 zwdlmtu(:,:) = 1._wp 164 zwdlmtv(:,:) = 1._wp 165 166 ! Horizontal Flux in u and v direction 167 DO jk = 1, jpkm1 168 DO jj = 1, jpj 169 DO ji = 1, jpi 170 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 171 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 172 END DO 173 END DO 174 END DO 175 176 zflxu(:,:) = zflxu(:,:) * e2u(:,:) 177 zflxv(:,:) = zflxv(:,:) * e1v(:,:) 178 179 wdmask(:,:) = 1 180 DO jj = 2, jpj 181 DO ji = 2, jpi 182 183 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE ! we don't care about land cells 184 IF( ht_0(ji,jj) - ssh_ref > zdepwd ) CYCLE ! and cells which are unlikely to dry 185 186 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & 187 & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji, jj-1), 0._wp) 188 zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj), 0._wp) + & 189 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 190 191 zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 192 IF(zdep2 .le. 0._wp) THEN !add more safty, but not necessary 193 sshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 194 IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp 195 IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 196 IF(zflxv(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = 0._wp 197 IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 198 wdmask(ji,jj) = 0._wp 199 END IF 200 ENDDO 201 END DO 202 203 204 ! HPG limiter from jholt 205 wdramp(:,:) = min((ht_0(:,:) + sshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp) 206 !jth assume don't need a lbc_lnk here 207 DO jj = 1, jpjm1 208 DO ji = 1, jpim1 209 wdrampu(ji,jj) = min(wdramp(ji,jj),wdramp(ji+1,jj)) 210 wdrampv(ji,jj) = min(wdramp(ji,jj),wdramp(ji,jj+1)) 211 END DO 212 END DO 213 ! end HPG limiter 214 181 215 182 216 183 217 !! start limiter iterations 184 185 186 187 188 189 190 191 192 218 DO jk1 = 1, nn_wdit + 1 219 220 221 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 222 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 223 jflag = 0 ! flag indicating if any further iterations are needed 224 225 DO jj = 2, jpj 226 DO ji = 2, jpi 193 227 194 195 IF( ht_wd(ji,jj) > zdepwd ) CYCLE228 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE 229 IF( ht_0(ji,jj) > zdepwd ) CYCLE 196 230 197 ztmp = e1e2t(ji,jj) 198 199 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) + & 200 & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp) 201 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + & 202 & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp) 203 204 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 205 zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 206 207 IF( zdep1 > zdep2 ) THEN 208 wdmask(ji, jj) = 0 209 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 210 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 211 ! flag if the limiter has been used but stop flagging if the only 212 ! changes have zeroed the coefficient since further iterations will 213 ! not change anything 214 IF( zcoef > 0._wp ) THEN 215 jflag = 1 216 ELSE 217 zcoef = 0._wp 218 ENDIF 219 IF(jk1 > nn_wdit) zcoef = 0._wp 220 IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef 221 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 222 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 223 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 224 END IF 225 END DO ! ji loop 226 END DO ! jj loop 227 228 CALL lbc_lnk( zwdlmtu, 'U', 1. ) 229 CALL lbc_lnk( zwdlmtv, 'V', 1. ) 230 231 IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain 232 233 IF(jflag == 0) EXIT 234 235 END DO ! jk1 loop 236 237 DO jk = 1, jpkm1 238 un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 239 vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 240 END DO 241 242 CALL lbc_lnk( un, 'U', -1. ) 243 CALL lbc_lnk( vn, 'V', -1. ) 244 ! 245 un_b(:,:) = un_b(:,:) * zwdlmtu(:, :) 246 vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :) 247 CALL lbc_lnk( un_b, 'U', -1. ) 248 CALL lbc_lnk( vn_b, 'V', -1. ) 249 250 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 251 252 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 253 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 231 ztmp = e1e2t(ji,jj) 232 233 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) + & 234 & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp) 235 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + & 236 & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp) 237 238 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 239 zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 240 241 IF( zdep1 > zdep2 ) THEN 242 wdmask(ji, jj) = 0 243 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 244 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 245 ! flag if the limiter has been used but stop flagging if the only 246 ! changes have zeroed the coefficient since further iterations will 247 ! not change anything 248 IF( zcoef > 0._wp ) THEN 249 jflag = 1 250 ELSE 251 zcoef = 0._wp 252 ENDIF 253 IF(jk1 > nn_wdit) zcoef = 0._wp 254 IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef 255 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 256 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 257 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 258 END IF 259 END DO ! ji loop 260 END DO ! jj loop 261 262 CALL lbc_lnk( zwdlmtu, 'U', 1. ) 263 CALL lbc_lnk( zwdlmtv, 'V', 1. ) 264 265 IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain 266 267 IF(jflag == 0) EXIT 268 269 END DO ! jk1 loop 270 271 DO jk = 1, jpkm1 272 un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 273 vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 274 END DO 275 276 CALL lbc_lnk( un, 'U', -1. ) 277 CALL lbc_lnk( vn, 'V', -1. ) 254 278 ! 255 ! 256 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 257 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 258 ! 259 ENDIF 279 un_b(:,:) = un_b(:,:) * zwdlmtu(:, :) 280 vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :) 281 CALL lbc_lnk( un_b, 'U', -1. ) 282 CALL lbc_lnk( vn_b, 'V', -1. ) 283 284 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 285 286 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 287 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 288 ! 289 ! 290 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 291 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 292 ! 260 293 ! 261 294 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') … … 291 324 IF( nn_timing == 1 ) CALL timing_start('wad_lmt_bt') 292 325 293 IF(ln_wd) THEN 294 295 CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 296 CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv) 297 ! 298 299 !IF(lwp) WRITE(numout,*) 300 !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting' 301 302 jflag = 0 303 zdepwd = 50._wp !maximum depth that ocean cells can have W/D processes 304 305 z2dt = rdtbt 306 307 zflxp(:,:) = 0._wp 308 zflxn(:,:) = 0._wp 309 310 zwdlmtu(:,:) = 1._wp 311 zwdlmtv(:,:) = 1._wp 326 327 CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 328 CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv) 329 ! 330 331 !IF(lwp) WRITE(numout,*) 332 !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting' 333 334 jflag = 0 335 zdepwd = 50._wp !maximum depth that ocean cells can have W/D processes 336 337 z2dt = rdtbt 338 339 zflxp(:,:) = 0._wp 340 zflxn(:,:) = 0._wp 341 342 zwdlmtu(:,:) = 1._wp 343 zwdlmtv(:,:) = 1._wp 312 344 313 345 ! Horizontal Flux in u and v direction 314 346 315 316 317 318 319 IF( ht_wd(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry320 321 322 323 324 325 326 zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1327 328 sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)329 330 331 332 333 334 335 347 DO jj = 2, jpj 348 DO ji = 2, jpi 349 350 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE ! we don't care about land cells 351 IF( ht_0(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry 352 353 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & 354 & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji, jj-1), 0._wp) 355 zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj), 0._wp) + & 356 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 357 358 zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 359 IF(zdep2 .le. 0._wp) THEN !add more safety, but not necessary 360 sshn_e(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 361 IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp 362 IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 363 IF(zflxv(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = 0._wp 364 IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 365 END IF 366 ENDDO 367 END DO 336 368 337 369 338 339 340 341 342 343 344 345 346 347 370 !! start limiter iterations 371 DO jk1 = 1, nn_wdit + 1 372 373 374 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 375 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 376 jflag = 0 ! flag indicating if any further iterations are needed 377 378 DO jj = 2, jpj 379 DO ji = 2, jpi 348 380 349 350 IF( ht_wd(ji,jj) > zdepwd ) CYCLE381 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE 382 IF( ht_0(ji,jj) > zdepwd ) CYCLE 351 383 352 ztmp = e1e2t(ji,jj) 353 354 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) + & 355 & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp) 356 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + & 357 & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp) 358 359 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 360 zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 361 362 IF(zdep1 > zdep2) THEN 363 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 364 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 365 ! flag if the limiter has been used but stop flagging if the only 366 ! changes have zeroed the coefficient since further iterations will 367 ! not change anything 368 IF( zcoef > 0._wp ) THEN 369 jflag = 1 370 ELSE 371 zcoef = 0._wp 372 ENDIF 373 IF(jk1 > nn_wdit) zcoef = 0._wp 374 IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef 375 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 376 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 377 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 378 END IF 379 END DO ! ji loop 380 END DO ! jj loop 381 382 CALL lbc_lnk( zwdlmtu, 'U', 1. ) 383 CALL lbc_lnk( zwdlmtv, 'V', 1. ) 384 385 IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain 386 387 IF(jflag == 0) EXIT 388 389 END DO ! jk1 loop 390 391 zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 392 zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 393 394 CALL lbc_lnk( zflxu, 'U', -1. ) 395 CALL lbc_lnk( zflxv, 'V', -1. ) 396 397 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 398 399 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 400 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 401 ! 402 ! 403 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 404 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 405 ! 406 END IF 384 ztmp = e1e2t(ji,jj) 385 386 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) + & 387 & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp) 388 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + & 389 & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp) 390 391 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 392 zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 393 394 IF(zdep1 > zdep2) THEN 395 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 396 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 397 ! flag if the limiter has been used but stop flagging if the only 398 ! changes have zeroed the coefficient since further iterations will 399 ! not change anything 400 IF( zcoef > 0._wp ) THEN 401 jflag = 1 402 ELSE 403 zcoef = 0._wp 404 ENDIF 405 IF(jk1 > nn_wdit) zcoef = 0._wp 406 IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef 407 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 408 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 409 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 410 END IF 411 END DO ! ji loop 412 END DO ! jj loop 413 414 CALL lbc_lnk( zwdlmtu, 'U', 1. ) 415 CALL lbc_lnk( zwdlmtv, 'V', 1. ) 416 417 IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain 418 419 IF(jflag == 0) EXIT 420 421 END DO ! jk1 loop 422 423 zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 424 zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 425 426 CALL lbc_lnk( zflxu, 'U', -1. ) 427 CALL lbc_lnk( zflxv, 'V', -1. ) 428 429 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 430 431 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 432 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 433 ! 434 ! 435 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 436 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 437 ! 407 438 408 439 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r6140 r8992 48 48 LOGICAL, PUBLIC :: ln_diaobs !: Logical switch for the obs operator 49 49 LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs 50 50 LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres 51 LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres 52 LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres 53 LOGICAL :: ln_sic_fp_indegs !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 54 55 REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres) 56 REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres) 57 REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres) 58 REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres) 59 REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres) 60 REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres) 61 REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres) 62 REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres) 63 51 64 INTEGER :: nn_1dint !: Vertical interpolation method 52 INTEGER :: nn_2dint !: Horizontal interpolation method 65 INTEGER :: nn_2dint !: Default horizontal interpolation method 66 INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method 67 INTEGER :: nn_2dint_sst !: SST horizontal interpolation method 68 INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method 69 INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method 53 70 INTEGER, DIMENSION(imaxavtypes) :: & 54 71 & nn_profdavtypes !: Profile data types representing a daily average … … 61 78 & nextrprof, & !: Number of profile extra variables 62 79 & nextrsurf !: Number of surface extra variables 63 INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: sstbias_type !SST bias type 80 INTEGER, DIMENSION(:), ALLOCATABLE :: & 81 & n2dintsurf !: Interpolation option for surface variables 82 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 83 & zavglamscl, & !: E/W diameter of averaging footprint for surface variables 84 & zavgphiscl !: N/S diameter of averaging footprint for surface variables 85 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 86 & lfpindegs, & !: T=> surface obs footprint size specified in degrees, F=> in metres 87 & llnightav !: Logical for calculating night-time averages 88 64 89 TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 65 90 & surfdata, & !: Initial surface data … … 116 141 & cn_profbfiles, & ! T/S profile input filenames 117 142 & cn_sstfbfiles, & ! Sea surface temperature input filenames 143 & cn_sssfbfiles, & ! Sea surface salinity input filenames 118 144 & cn_slafbfiles, & ! Sea level anomaly input filenames 119 145 & cn_sicfbfiles, & ! Seaice concentration input filenames 120 146 & cn_velfbfiles, & ! Velocity profile input filenames 121 & cn_sstbias _files ! SST bias input filenames147 & cn_sstbiasfiles ! SST bias input filenames 122 148 CHARACTER(LEN=128) :: & 123 149 & cn_altbiasfile ! Altimeter bias input filename … … 130 156 LOGICAL :: ln_sla ! Logical switch for sea level anomalies 131 157 LOGICAL :: ln_sst ! Logical switch for sea surface temperature 158 LOGICAL :: ln_sss ! Logical switch for sea surface salinity 132 159 LOGICAL :: ln_sic ! Logical switch for sea ice concentration 133 160 LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs 134 161 LOGICAL :: ln_nea ! Logical switch to remove obs near land 135 162 LOGICAL :: ln_altbias ! Logical switch for altimeter bias 136 LOGICAL :: ln_sstbias !:Logical switch for bias corection of SST163 LOGICAL :: ln_sstbias ! Logical switch for bias corection of SST 137 164 LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files 138 165 LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs 166 LOGICAL :: ln_bound_reject ! Logical to remove obs near boundaries in LAMs. 139 167 LOGICAL :: llvar1 ! Logical for profile variable 1 140 168 LOGICAL :: llvar2 ! Logical for profile variable 1 141 LOGICAL :: llnightav ! Logical for calculating night-time averages142 169 LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files 143 170 … … 155 182 156 183 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & 157 & ln_sst, ln_sic, ln_vel3d, & 158 & ln_altbias, ln_nea, ln_grid_global, & 159 & ln_grid_search_lookup, & 160 & ln_ignmis, ln_s_at_t, ln_sstnight, & 184 & ln_sst, ln_sic, ln_sss, ln_vel3d, & 185 & ln_altbias, ln_sstbias, ln_nea, & 186 & ln_grid_global, ln_grid_search_lookup, & 187 & ln_ignmis, ln_s_at_t, ln_bound_reject, & 188 & ln_sstnight, & 189 & ln_sla_fp_indegs, ln_sst_fp_indegs, & 190 & ln_sss_fp_indegs, ln_sic_fp_indegs, & 161 191 & cn_profbfiles, cn_slafbfiles, & 162 192 & cn_sstfbfiles, cn_sicfbfiles, & 163 & cn_velfbfiles, cn_altbiasfile, & 193 & cn_velfbfiles, cn_sssfbfiles, & 194 & cn_sstbiasfiles, cn_altbiasfile, & 164 195 & cn_gridsearchfile, rn_gridsearchres, & 165 & rn_dobsini, rn_dobsend, nn_1dint, nn_2dint, & 196 & rn_dobsini, rn_dobsend, & 197 & rn_sla_avglamscl, rn_sla_avgphiscl, & 198 & rn_sst_avglamscl, rn_sst_avgphiscl, & 199 & rn_sss_avglamscl, rn_sss_avgphiscl, & 200 & rn_sic_avglamscl, rn_sic_avgphiscl, & 201 & nn_1dint, nn_2dint, & 202 & nn_2dint_sla, nn_2dint_sst, & 203 & nn_2dint_sss, nn_2dint_sic, & 166 204 & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & 167 & nn_profdavtypes , ln_sstbias, cn_sstbias_files205 & nn_profdavtypes 168 206 169 207 INTEGER :: jnumsstbias … … 178 216 ! Read namelist parameters 179 217 !----------------------------------------------------------------------- 180 181 !Initalise all values in namelist arrays182 ALLOCATE(sstbias_type(jpmaxnfiles))183 218 ! Some namelist arrays need initialising 184 219 cn_profbfiles(:) = '' … … 187 222 cn_sicfbfiles(:) = '' 188 223 cn_velfbfiles(:) = '' 189 cn_sstbias_files(:) = '' 224 cn_sssfbfiles(:) = '' 225 cn_sstbiasfiles(:) = '' 190 226 nn_profdavtypes(:) = -1 191 227 … … 208 244 RETURN 209 245 ENDIF 210 211 !----------------------------------------------------------------------- 212 ! Set up list of observation types to be used 213 ! and the files associated with each type 214 !----------------------------------------------------------------------- 215 216 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 217 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) ) 218 219 IF (ln_sstbias) THEN 220 lmask(:) = .FALSE. 221 WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE. 222 jnumsstbias = COUNT(lmask) 223 lmask(:) = .FALSE. 224 ENDIF 225 226 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 227 IF(lwp) WRITE(numout,cform_war) 228 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 229 & ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 230 & ' are set to .FALSE. so turning off calls to dia_obs' 231 nwarn = nwarn + 1 232 ln_diaobs = .FALSE. 233 RETURN 234 ENDIF 235 236 IF ( nproftypes > 0 ) THEN 237 238 ALLOCATE( cobstypesprof(nproftypes) ) 239 ALLOCATE( ifilesprof(nproftypes) ) 240 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 241 242 jtype = 0 243 IF (ln_t3d .OR. ln_s3d) THEN 244 jtype = jtype + 1 245 clproffiles(jtype,:) = cn_profbfiles(:) 246 cobstypesprof(jtype) = 'prof ' 247 ifilesprof(jtype) = 0 248 DO jfile = 1, jpmaxnfiles 249 IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 250 ifilesprof(jtype) = ifilesprof(jtype) + 1 251 END DO 252 ENDIF 253 IF (ln_vel3d) THEN 254 jtype = jtype + 1 255 clproffiles(jtype,:) = cn_velfbfiles(:) 256 cobstypesprof(jtype) = 'vel ' 257 ifilesprof(jtype) = 0 258 DO jfile = 1, jpmaxnfiles 259 IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 260 ifilesprof(jtype) = ifilesprof(jtype) + 1 261 END DO 262 ENDIF 263 264 ENDIF 265 266 IF ( nsurftypes > 0 ) THEN 267 268 ALLOCATE( cobstypessurf(nsurftypes) ) 269 ALLOCATE( ifilessurf(nsurftypes) ) 270 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 271 272 jtype = 0 273 IF (ln_sla) THEN 274 jtype = jtype + 1 275 clsurffiles(jtype,:) = cn_slafbfiles(:) 276 cobstypessurf(jtype) = 'sla ' 277 ifilessurf(jtype) = 0 278 DO jfile = 1, jpmaxnfiles 279 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 280 ifilessurf(jtype) = ifilessurf(jtype) + 1 281 END DO 282 ENDIF 283 IF (ln_sst) THEN 284 jtype = jtype + 1 285 clsurffiles(jtype,:) = cn_sstfbfiles(:) 286 cobstypessurf(jtype) = 'sst ' 287 ifilessurf(jtype) = 0 288 DO jfile = 1, jpmaxnfiles 289 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 290 ifilessurf(jtype) = ifilessurf(jtype) + 1 291 END DO 292 ENDIF 293 #if defined key_lim2 || defined key_lim3 294 IF (ln_sic) THEN 295 jtype = jtype + 1 296 clsurffiles(jtype,:) = cn_sicfbfiles(:) 297 cobstypessurf(jtype) = 'sic ' 298 ifilessurf(jtype) = 0 299 DO jfile = 1, jpmaxnfiles 300 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 301 ifilessurf(jtype) = ifilessurf(jtype) + 1 302 END DO 303 ENDIF 304 #endif 305 306 ENDIF 307 308 !Write namelist settings to stdout 246 309 247 IF(lwp) THEN 310 248 WRITE(numout,*) … … 318 256 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic 319 257 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 320 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ',ln_grid_global321 WRITE(numout,*) ' Logical switch for SST bias correction ln_sstbias = ', ln_sstbias322 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup258 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 259 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global 260 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 323 261 IF (ln_grid_search_lookup) & 324 262 WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile … … 328 266 WRITE(numout,*) ' Type of horizontal interpolation method nn_2dint = ', nn_2dint 329 267 WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea 268 WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject 330 269 WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc 331 270 WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr 332 271 WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff 333 272 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 273 WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias 334 274 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 335 275 WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes 336 276 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 337 WRITE(numout,*) ' Number of profile obs types: ',nproftypes 338 339 IF ( nproftypes > 0 ) THEN 340 DO jtype = 1, nproftypes 341 DO jfile = 1, ifilesprof(jtype) 342 WRITE(numout,'(1X,2A)') ' '//cobstypesprof(jtype)//' input observation file names = ', & 343 TRIM(clproffiles(jtype,jfile)) 344 END DO 345 END DO 277 ENDIF 278 !----------------------------------------------------------------------- 279 ! Set up list of observation types to be used 280 ! and the files associated with each type 281 !----------------------------------------------------------------------- 282 283 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 284 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss /) ) 285 286 IF (ln_sstbias) THEN 287 lmask(:) = .FALSE. 288 WHERE (cn_sstbiasfiles(:) /= '') lmask(:) = .TRUE. 289 jnumsstbias = COUNT(lmask) 290 lmask(:) = .FALSE. 291 ENDIF 292 293 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 294 IF(lwp) WRITE(numout,cform_war) 295 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 296 & ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 297 & ' are set to .FALSE. so turning off calls to dia_obs' 298 nwarn = nwarn + 1 299 ln_diaobs = .FALSE. 300 RETURN 301 ENDIF 302 303 IF ( nproftypes > 0 ) THEN 304 305 ALLOCATE( cobstypesprof(nproftypes) ) 306 ALLOCATE( ifilesprof(nproftypes) ) 307 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 308 309 jtype = 0 310 IF (ln_t3d .OR. ln_s3d) THEN 311 jtype = jtype + 1 312 CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof ', & 313 & cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 346 314 ENDIF 347 348 WRITE(numout,*)' Number of surface obs types: ',nsurftypes 349 IF ( nsurftypes > 0 ) THEN 350 DO jtype = 1, nsurftypes 351 DO jfile = 1, ifilessurf(jtype) 352 WRITE(numout,'(1X,2A)') ' '//cobstypessurf(jtype)//' input observation file names = ', & 353 TRIM(clsurffiles(jtype,jfile)) 354 END DO 355 END DO 315 IF (ln_vel3d) THEN 316 jtype = jtype + 1 317 CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel ', & 318 & cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 356 319 ENDIF 357 WRITE(numout,*) '~~~~~~~~~~~~' 358 359 ENDIF 320 321 ENDIF 322 323 IF ( nsurftypes > 0 ) THEN 324 325 ALLOCATE( cobstypessurf(nsurftypes) ) 326 ALLOCATE( ifilessurf(nsurftypes) ) 327 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 328 ALLOCATE(n2dintsurf(nsurftypes)) 329 ALLOCATE(zavglamscl(nsurftypes)) 330 ALLOCATE(zavgphiscl(nsurftypes)) 331 ALLOCATE(lfpindegs(nsurftypes)) 332 ALLOCATE(llnightav(nsurftypes)) 333 334 jtype = 0 335 IF (ln_sla) THEN 336 jtype = jtype + 1 337 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla ', & 338 & cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles ) 339 CALL obs_setinterpopts( nsurftypes, jtype, 'sla ', & 340 & nn_2dint, nn_2dint_sla, & 341 & rn_sla_avglamscl, rn_sla_avgphiscl, & 342 & ln_sla_fp_indegs, .FALSE., & 343 & n2dintsurf, zavglamscl, zavgphiscl, & 344 & lfpindegs, llnightav ) 345 ENDIF 346 IF (ln_sst) THEN 347 jtype = jtype + 1 348 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst ', & 349 & cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 350 CALL obs_setinterpopts( nsurftypes, jtype, 'sst ', & 351 & nn_2dint, nn_2dint_sst, & 352 & rn_sst_avglamscl, rn_sst_avgphiscl, & 353 & ln_sst_fp_indegs, ln_sstnight, & 354 & n2dintsurf, zavglamscl, zavgphiscl, & 355 & lfpindegs, llnightav ) 356 ENDIF 357 #if defined key_lim2 || defined key_lim3 || defined key_cice 358 IF (ln_sic) THEN 359 jtype = jtype + 1 360 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic ', & 361 & cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 362 CALL obs_setinterpopts( nsurftypes, jtype, 'sic ', & 363 & nn_2dint, nn_2dint_sic, & 364 & rn_sic_avglamscl, rn_sic_avgphiscl, & 365 & ln_sic_fp_indegs, .FALSE., & 366 & n2dintsurf, zavglamscl, zavgphiscl, & 367 & lfpindegs, llnightav ) 368 ENDIF 369 #endif 370 IF (ln_sss) THEN 371 jtype = jtype + 1 372 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss ', & 373 & cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 374 CALL obs_setinterpopts( nsurftypes, jtype, 'sss ', & 375 & nn_2dint, nn_2dint_sss, & 376 & rn_sss_avglamscl, rn_sss_avgphiscl, & 377 & ln_sss_fp_indegs, .FALSE., & 378 & n2dintsurf, zavglamscl, zavgphiscl, & 379 & lfpindegs, llnightav ) 380 ENDIF 381 382 ENDIF 383 384 360 385 361 386 !----------------------------------------------------------------------- … … 377 402 ENDIF 378 403 379 IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4) ) THEN404 IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN 380 405 CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 381 406 & ' is not available') … … 442 467 & jpi, jpj, jpk, & 443 468 & zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, & 444 & ln_nea, kdailyavtypes = nn_profdavtypes ) 469 & ln_nea, ln_bound_reject, & 470 & kdailyavtypes = nn_profdavtypes ) 445 471 446 472 END DO … … 461 487 nvarssurf(jtype) = 1 462 488 nextrsurf(jtype) = 0 463 llnightav = .FALSE.489 llnightav(jtype) = .FALSE. 464 490 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 465 IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight491 IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav(jtype) = ln_sstnight 466 492 467 493 !Read in surface obs types … … 469 495 & clsurffiles(jtype,1:ifilessurf(jtype)), & 470 496 & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 471 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 472 473 474 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 497 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 498 499 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 475 500 476 501 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 477 CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 478 IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 502 CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 503 IF ( ln_altbias ) & 504 & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 479 505 ENDIF 506 480 507 IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 481 !Read in bias field and correct SST. 482 IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 483 " but no bias"// & 484 " files to read in") 485 CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, & 486 jnumsstbias, cn_sstbias_files(1:jnumsstbias) ) 508 jnumsstbias = 0 509 DO jfile = 1, jpmaxnfiles 510 IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 511 & jnumsstbias = jnumsstbias + 1 512 END DO 513 IF ( jnumsstbias == 0 ) THEN 514 CALL ctl_stop("ln_sstbias set but no bias files to read in") 515 ENDIF 516 517 CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 518 & jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 519 487 520 ENDIF 488 521 END DO … … 545 578 & frld 546 579 #endif 580 #if defined key_cice 581 USE sbc_oce, ONLY : fr_i ! ice fraction 582 #endif 583 547 584 IMPLICIT NONE 548 585 … … 561 598 & zprofmask2 ! Mask associated with zprofvar2 562 599 REAL(wp), POINTER, DIMENSION(:,:) :: & 563 & zsurfvar ! Model values equivalent to surface ob. 600 & zsurfvar, & ! Model values equivalent to surface ob. 601 & zsurfmask ! Mask associated with surface variable 564 602 REAL(wp), POINTER, DIMENSION(:,:) :: & 565 603 & zglam1, & ! Model longitudes for prof variable 1 … … 567 605 & zgphi1, & ! Model latitudes for prof variable 1 568 606 & zgphi2 ! Model latitudes for prof variable 2 569 #if ! defined key_lim2 && ! defined key_lim3570 REAL(wp), POINTER, DIMENSION(:,:) :: frld571 #endif572 LOGICAL :: llnightav ! Logical for calculating night-time average573 607 574 608 !Allocate local work arrays … … 578 612 CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 579 613 CALL wrk_alloc( jpi, jpj, zsurfvar ) 614 CALL wrk_alloc( jpi, jpj, zsurfmask ) 580 615 CALL wrk_alloc( jpi, jpj, zglam1 ) 581 616 CALL wrk_alloc( jpi, jpj, zglam2 ) 582 617 CALL wrk_alloc( jpi, jpj, zgphi1 ) 583 618 CALL wrk_alloc( jpi, jpj, zgphi2 ) 584 #if ! defined key_lim2 && ! defined key_lim3585 CALL wrk_alloc(jpi,jpj,frld)586 #endif587 619 588 620 IF(lwp) THEN … … 594 626 idaystp = NINT( rday / rdt ) 595 627 596 !-----------------------------------------------------------------------597 ! No LIM => frld == 0.0_wp598 !-----------------------------------------------------------------------599 #if ! defined key_lim2 && ! defined key_lim3600 frld(:,:) = 0.0_wp601 #endif602 628 !----------------------------------------------------------------------- 603 629 ! Call the profile and surface observation operators … … 627 653 zgphi1(:,:) = gphiu(:,:) 628 654 zgphi2(:,:) = gphiv(:,:) 655 CASE DEFAULT 656 CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 629 657 END SELECT 630 658 631 IF( ln_zco .OR. ln_zps ) THEN 632 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 633 & nit000, idaystp, & 634 & zprofvar1, zprofvar2, & 635 & gdept_1d, zprofmask1, zprofmask2, & 636 & zglam1, zglam2, zgphi1, zgphi2, & 637 & nn_1dint, nn_2dint, & 638 & kdailyavtypes = nn_profdavtypes ) 639 ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') THEN 640 !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION 641 CALL obs_pro_sco_opt( profdataqc(jtype), & 642 & kstp, jpi, jpj, jpk, nit000, idaystp, & 643 & zprofvar1, zprofvar2, & 644 & gdept_n(:,:,:), gdepw_n(:,:,:), & 645 & tmask, nn_1dint, nn_2dint, & 646 & kdailyavtypes = nn_profdavtypes ) 647 ELSE 648 CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// & 649 'yet working for velocity data (turn off velocity observations') 650 ENDIF 659 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 660 & nit000, idaystp, & 661 & zprofvar1, zprofvar2, & 662 & gdept_n(:,:,:), gdepw_n(:,:,:), & 663 & zprofmask1, zprofmask2, & 664 & zglam1, zglam2, zgphi1, zgphi2, & 665 & nn_1dint, nn_2dint, & 666 & kdailyavtypes = nn_profdavtypes ) 651 667 652 668 END DO … … 657 673 658 674 DO jtype = 1, nsurftypes 675 676 !Defaults which might be changed 677 zsurfmask(:,:) = tmask(:,:,1) 659 678 660 679 SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 661 680 CASE('sst') 662 681 zsurfvar(:,:) = tsn(:,:,1,jp_tem) 663 llnightav = ln_sstnight664 682 CASE('sla') 665 683 zsurfvar(:,:) = sshn(:,:) 666 llnightav = .FALSE.667 #if defined key_lim2 || defined key_lim3 684 CASE('sss') 685 zsurfvar(:,:) = tsn(:,:,1,jp_sal) 668 686 CASE('sic') 669 687 IF ( kstp == 0 ) THEN … … 678 696 CYCLE 679 697 ELSE 698 #if defined key_cice 699 zsurfvar(:,:) = fr_i(:,:) 700 #elif defined key_lim2 || defined key_lim3 680 701 zsurfvar(:,:) = 1._wp - frld(:,:) 702 #else 703 CALL ctl_stop( ' Trying to run sea-ice observation operator', & 704 & ' but no sea-ice model appears to have been defined' ) 705 #endif 681 706 ENDIF 682 707 683 llnightav = .FALSE.684 #endif685 708 END SELECT 686 709 687 710 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 688 & nit000, idaystp, zsurfvar, tmask(:,:,1), & 689 & nn_2dint, llnightav ) 711 & nit000, idaystp, zsurfvar, zsurfmask, & 712 & n2dintsurf(jtype), llnightav(jtype), & 713 & zavglamscl(jtype), zavgphiscl(jtype), & 714 & lfpindegs(jtype) ) 690 715 691 716 END DO … … 698 723 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 699 724 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 725 CALL wrk_dealloc( jpi, jpj, zsurfmask ) 700 726 CALL wrk_dealloc( jpi, jpj, zglam1 ) 701 727 CALL wrk_dealloc( jpi, jpj, zglam2 ) 702 728 CALL wrk_dealloc( jpi, jpj, zgphi1 ) 703 729 CALL wrk_dealloc( jpi, jpj, zgphi2 ) 704 #if ! defined key_lim2 && ! defined key_lim3705 CALL wrk_dealloc(jpi,jpj,frld)706 #endif707 730 708 731 END SUBROUTINE dia_obs … … 821 844 822 845 IF ( nsurftypes > 0 ) & 823 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 846 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 847 & n2dintsurf, zavglamscl, zavgphiscl, lfpindegs, llnightav ) 824 848 825 849 END SUBROUTINE dia_obs_dealloc … … 970 994 END SUBROUTINE fin_date 971 995 996 SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, & 997 & cfilestype, ifiles, cobstypes, cfiles ) 998 999 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 1000 INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 1001 INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs 1002 INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 1003 & ifiles ! Out appended number of files for this type 1004 1005 CHARACTER(len=6), INTENT(IN) :: ctypein 1006 CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: & 1007 & cfilestype ! In list of files for this obs type 1008 CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: & 1009 & cobstypes ! Out appended list of obs types 1010 CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: & 1011 & cfiles ! Out appended list of files for all types 1012 1013 !Local variables 1014 INTEGER :: jfile 1015 1016 cfiles(jtype,:) = cfilestype(:) 1017 cobstypes(jtype) = ctypein 1018 ifiles(jtype) = 0 1019 DO jfile = 1, jpmaxnfiles 1020 IF ( trim(cfiles(jtype,jfile)) /= '' ) & 1021 ifiles(jtype) = ifiles(jtype) + 1 1022 END DO 1023 1024 IF ( ifiles(jtype) == 0 ) THEN 1025 CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)// & 1026 & ' set to true but no files available to read' ) 1027 ENDIF 1028 1029 IF(lwp) THEN 1030 WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:' 1031 DO jfile = 1, ifiles(jtype) 1032 WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile)) 1033 END DO 1034 ENDIF 1035 1036 END SUBROUTINE obs_settypefiles 1037 1038 SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, & 1039 & n2dint_default, n2dint_type, & 1040 & zavglamscl_type, zavgphiscl_type, & 1041 & lfp_indegs_type, lavnight_type, & 1042 & n2dint, zavglamscl, zavgphiscl, & 1043 & lfpindegs, lavnight ) 1044 1045 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 1046 INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs 1047 INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type 1048 INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type 1049 REAL(wp), INTENT(IN) :: & 1050 & zavglamscl_type, & !E/W diameter of obs footprint for this type 1051 & zavgphiscl_type !N/S diameter of obs footprint for this type 1052 LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres 1053 LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average 1054 CHARACTER(len=6), INTENT(IN) :: ctypein 1055 1056 INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 1057 & n2dint 1058 REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 1059 & zavglamscl, zavgphiscl 1060 LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 1061 & lfpindegs, lavnight 1062 1063 lavnight(jtype) = lavnight_type 1064 1065 IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN 1066 n2dint(jtype) = n2dint_type 1067 ELSE 1068 n2dint(jtype) = n2dint_default 1069 ENDIF 1070 1071 ! For averaging observation footprints set options for size of footprint 1072 IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 1073 IF ( zavglamscl_type > 0._wp ) THEN 1074 zavglamscl(jtype) = zavglamscl_type 1075 ELSE 1076 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 1077 'scale (zavglamscl) for observation type '//TRIM(ctypein) ) 1078 ENDIF 1079 1080 IF ( zavgphiscl_type > 0._wp ) THEN 1081 zavgphiscl(jtype) = zavgphiscl_type 1082 ELSE 1083 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 1084 'scale (zavgphiscl) for observation type '//TRIM(ctypein) ) 1085 ENDIF 1086 1087 lfpindegs(jtype) = lfp_indegs_type 1088 1089 ENDIF 1090 1091 ! Write out info 1092 IF(lwp) THEN 1093 IF ( n2dint(jtype) <= 4 ) THEN 1094 WRITE(numout,*) ' '//TRIM(ctypein)// & 1095 & ' model counterparts will be interpolated horizontally' 1096 ELSE IF ( n2dint(jtype) <= 6 ) THEN 1097 WRITE(numout,*) ' '//TRIM(ctypein)// & 1098 & ' model counterparts will be averaged horizontally' 1099 WRITE(numout,*) ' '//' with E/W scale: ',zavglamscl(jtype) 1100 WRITE(numout,*) ' '//' with N/S scale: ',zavgphiscl(jtype) 1101 IF ( lfpindegs(jtype) ) THEN 1102 WRITE(numout,*) ' '//' (in degrees)' 1103 ELSE 1104 WRITE(numout,*) ' '//' (in metres)' 1105 ENDIF 1106 ENDIF 1107 ENDIF 1108 1109 END SUBROUTINE obs_setinterpopts 1110 972 1111 END MODULE diaobs -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_mpp.F90
r6140 r8992 142 142 !! 143 143 144 iobsp =kobsp144 iobsp(:)=kobsp(:) 145 145 146 146 WHERE( iobsp(:) == -1 ) … … 148 148 END WHERE 149 149 150 iobsp =-1*iobsp150 iobsp(:)=-1*iobsp(:) 151 151 152 152 CALL obs_mpp_max_integer( iobsp, kno ) 153 153 154 kobsp =-1*iobsp154 kobsp(:)=-1*iobsp(:) 155 155 156 156 isum=0 -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r7646 r8992 9 9 !! obs_prof_opt : Compute the model counterpart of profile data 10 10 !! obs_surf_opt : Compute the model counterpart of surface data 11 !! obs_pro_sco_opt: Compute the model counterpart of temperature and12 !! salinity observations from profiles in generalised13 !! vertical coordinates14 11 !!---------------------------------------------------------------------- 15 12 … … 22 19 & obs_int_h2d, & 23 20 & obs_int_h2d_init 21 USE obs_averg_h2d, ONLY : & ! Horizontal averaging to the obs footprint 22 & obs_avg_h2d, & 23 & obs_avg_h2d_init, & 24 & obs_max_fpsize 24 25 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt 25 26 & obs_int_z1d, & 26 27 & obs_int_z1d_spl 27 USE obs_const, ONLY : &28 & obfillflt ! Fillvalue28 USE obs_const, ONLY : & ! Obs fill value 29 & obfillflt 29 30 USE dom_oce, ONLY : & 30 & glamt, glamu, glamv, & 31 & gphit, gphiu, gphiv, & 32 & gdept_n, gdept_0 33 USE lib_mpp, ONLY : & 31 & glamt, glamf, & 32 & gphit, gphif 33 USE lib_mpp, ONLY : & ! Warning and stopping routines 34 34 & ctl_warn, ctl_stop 35 USE sbcdcy, ONLY : & ! For calculation of where it is night-time 36 & sbc_dcy, nday_qsr 35 37 USE obs_grid, ONLY : & 36 38 & obs_level_search 37 USE sbcdcy, ONLY : & ! For calculation of where it is night-time38 & sbc_dcy, nday_qsr39 39 40 40 IMPLICIT NONE … … 44 44 45 45 PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs 46 & obs_pro_sco_opt, & ! Compute the model counterpart of profile observations47 46 & obs_surf_opt ! Compute the model counterpart of surface obs 48 47 … … 58 57 CONTAINS 59 58 59 60 60 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 61 61 & kit000, kdaystp, & 62 & pvar1, pvar2, pgdept, pmask1, pmask2, & 62 & pvar1, pvar2, pgdept, pgdepw, & 63 & pmask1, pmask2, & 63 64 & plam1, plam2, pphi1, pphi2, & 64 65 & k1dint, k2dint, kdailyavtypes ) … … 111 112 !! ! 07-03 (K. Mogensen) General handling of profiles 112 113 !! ! 15-02 (M. Martin) Combined routine for all profile types 114 !! ! 17-02 (M. Martin) Include generalised vertical coordinate changes 113 115 !!----------------------------------------------------------------------- 114 116 … … 140 142 & pphi1, & ! Model latitudes for variable 1 141 143 & pphi2 ! Model latitudes for variable 2 142 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 143 & pgdept ! Model array of depth levels 144 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 145 & pgdept, & ! Model array of depth T levels 146 & pgdepw ! Model array of depth W levels 144 147 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 145 148 & kdailyavtypes ! Types for daily averages … … 156 159 INTEGER :: iend 157 160 INTEGER :: iobs 161 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 162 INTEGER :: inum_obs 158 163 INTEGER, DIMENSION(imaxavtypes) :: & 159 164 & idailyavtypes … … 163 168 & igrdj1, & 164 169 & igrdj2 170 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 171 165 172 REAL(KIND=wp) :: zlam 166 173 REAL(KIND=wp) :: zphi … … 171 178 & zobsk, & 172 179 & zobs2k 173 REAL(KIND=wp), DIMENSION(2,2, kpk) :: &180 REAL(KIND=wp), DIMENSION(2,2,1) :: & 174 181 & zweig1, & 175 & zweig2 182 & zweig2, & 183 & zweig 176 184 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 177 185 & zmask1, & 178 186 & zmask2, & 179 & zint1, & 180 & zint2, & 181 & zinm1, & 182 & zinm2 187 & zint1, & 188 & zint2, & 189 & zinm1, & 190 & zinm2, & 191 & zgdept, & 192 & zgdepw 183 193 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 184 194 & zglam1, & … … 186 196 & zgphi1, & 187 197 & zgphi2 198 REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2 199 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 200 188 201 LOGICAL :: ld_dailyav 189 202 … … 266 279 & zmask1(2,2,kpk,ipro), & 267 280 & zmask2(2,2,kpk,ipro), & 268 & zint1(2,2,kpk,ipro), & 269 & zint2(2,2,kpk,ipro) & 281 & zint1(2,2,kpk,ipro), & 282 & zint2(2,2,kpk,ipro), & 283 & zgdept(2,2,kpk,ipro), & 284 & zgdepw(2,2,kpk,ipro) & 270 285 & ) 271 286 … … 290 305 END DO 291 306 307 ! Initialise depth arrays 308 zgdept(:,:,:,:) = 0.0 309 zgdepw(:,:,:,:) = 0.0 310 292 311 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 293 312 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) … … 300 319 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 ) 301 320 321 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 322 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 323 302 324 ! At the end of the day also get interpolated means 303 325 IF ( ld_dailyav .AND. idayend == 0 ) THEN … … 314 336 315 337 ENDIF 338 339 ! Return if no observations to process 340 ! Has to be done after comm commands to ensure processors 341 ! stay in sync 342 IF ( ipro == 0 ) RETURN 316 343 317 344 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro … … 339 366 zphi = prodatqc%rphi(jobs) 340 367 341 ! Horizontal weights and vertical mask342 368 ! Horizontal weights 369 ! Masked values are calculated later. 343 370 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 344 371 345 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &372 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 346 373 & zglam1(:,:,iobs), zgphi1(:,:,iobs), & 347 & zmask1(:,:, :,iobs), zweig1, zobsmask1 )374 & zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 348 375 349 376 ENDIF … … 351 378 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 352 379 353 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &380 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 354 381 & zglam2(:,:,iobs), zgphi2(:,:,iobs), & 355 & zmask2(:,:, :,iobs), zweig2, zobsmask2 )382 & zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 356 383 357 384 ENDIF … … 365 392 IF ( idayend == 0 ) THEN 366 393 ! Daily averaged data 367 CALL obs_int_h2d( kpk, kpk, & 368 & zweig1, zinm1(:,:,:,iobs), zobsk ) 369 370 ENDIF 371 372 ELSE 373 374 ! Point data 375 CALL obs_int_h2d( kpk, kpk, & 376 & zweig1, zint1(:,:,:,iobs), zobsk ) 377 378 ENDIF 379 380 !------------------------------------------------------------- 381 ! Compute vertical second-derivative of the interpolating 382 ! polynomial at obs points 383 !------------------------------------------------------------- 384 385 IF ( k1dint == 1 ) THEN 386 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 387 & pgdept, zobsmask1 ) 388 ENDIF 389 390 !----------------------------------------------------------------- 391 ! Vertical interpolation to the observation point 392 !----------------------------------------------------------------- 393 ista = prodatqc%npvsta(jobs,1) 394 iend = prodatqc%npvend(jobs,1) 395 CALL obs_int_z1d( kpk, & 396 & prodatqc%var(1)%mvk(ista:iend), & 397 & k1dint, iend - ista + 1, & 398 & prodatqc%var(1)%vdep(ista:iend), & 399 & zobsk, zobs2k, & 400 & prodatqc%var(1)%vmod(ista:iend), & 401 & pgdept, zobsmask1 ) 402 403 ENDIF 404 394 395 ! vertically interpolate all 4 corners 396 ista = prodatqc%npvsta(jobs,1) 397 iend = prodatqc%npvend(jobs,1) 398 inum_obs = iend - ista + 1 399 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 400 401 DO iin=1,2 402 DO ijn=1,2 403 404 IF ( k1dint == 1 ) THEN 405 CALL obs_int_z1d_spl( kpk, & 406 & zinm1(iin,ijn,:,iobs), & 407 & zobs2k, zgdept(iin,ijn,:,iobs), & 408 & zmask1(iin,ijn,:,iobs)) 409 ENDIF 410 411 CALL obs_level_search(kpk, & 412 & zgdept(iin,ijn,:,iobs), & 413 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 414 & iv_indic) 415 416 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 417 & prodatqc%var(1)%vdep(ista:iend), & 418 & zinm1(iin,ijn,:,iobs), & 419 & zobs2k, interp_corner(iin,ijn,:), & 420 & zgdept(iin,ijn,:,iobs), & 421 & zmask1(iin,ijn,:,iobs)) 422 423 ENDDO 424 ENDDO 425 426 ENDIF !idayend 427 428 ELSE 429 430 ! Point data 431 432 ! vertically interpolate all 4 corners 433 ista = prodatqc%npvsta(jobs,1) 434 iend = prodatqc%npvend(jobs,1) 435 inum_obs = iend - ista + 1 436 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 437 DO iin=1,2 438 DO ijn=1,2 439 440 IF ( k1dint == 1 ) THEN 441 CALL obs_int_z1d_spl( kpk, & 442 & zint1(iin,ijn,:,iobs),& 443 & zobs2k, zgdept(iin,ijn,:,iobs), & 444 & zmask1(iin,ijn,:,iobs)) 445 446 ENDIF 447 448 CALL obs_level_search(kpk, & 449 & zgdept(iin,ijn,:,iobs),& 450 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 451 & iv_indic) 452 453 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 454 & prodatqc%var(1)%vdep(ista:iend), & 455 & zint1(iin,ijn,:,iobs), & 456 & zobs2k,interp_corner(iin,ijn,:), & 457 & zgdept(iin,ijn,:,iobs), & 458 & zmask1(iin,ijn,:,iobs) ) 459 460 ENDDO 461 ENDDO 462 463 ENDIF 464 465 !------------------------------------------------------------- 466 ! Compute the horizontal interpolation for every profile level 467 !------------------------------------------------------------- 468 469 DO ikn=1,inum_obs 470 iend=ista+ikn-1 471 472 zweig(:,:,1) = 0._wp 473 474 ! This code forces the horizontal weights to be 475 ! zero IF the observation is below the bottom of the 476 ! corners of the interpolation nodes, Or if it is in 477 ! the mask. This is important for observations near 478 ! steep bathymetry 479 DO iin=1,2 480 DO ijn=1,2 481 482 depth_loop1: DO ik=kpk,2,-1 483 IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN 484 485 zweig(iin,ijn,1) = & 486 & zweig1(iin,ijn,1) * & 487 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 488 & - prodatqc%var(1)%vdep(iend)),0._wp) 489 490 EXIT depth_loop1 491 492 ENDIF 493 494 ENDDO depth_loop1 495 496 ENDDO 497 ENDDO 498 499 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 500 & prodatqc%var(1)%vmod(iend:iend) ) 501 502 ! Set QC flag for any observations found below the bottom 503 ! needed as the check here is more strict than that in obs_prep 504 IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4 505 506 ENDDO 507 508 DEALLOCATE(interp_corner,iv_indic) 509 510 ENDIF 511 512 ! For the second variable 405 513 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 406 514 … … 410 518 411 519 IF ( idayend == 0 ) THEN 412 413 520 ! Daily averaged data 414 CALL obs_int_h2d( kpk, kpk, & 415 & zweig2, zinm2(:,:,:,iobs), zobsk ) 416 417 ENDIF 418 419 ELSE 420 421 ! Point data 422 CALL obs_int_h2d( kpk, kpk, & 423 & zweig2, zint2(:,:,:,iobs), zobsk ) 424 425 ENDIF 426 427 428 !------------------------------------------------------------- 429 ! Compute vertical second-derivative of the interpolating 430 ! polynomial at obs points 431 !------------------------------------------------------------- 432 433 IF ( k1dint == 1 ) THEN 434 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 435 & pgdept, zobsmask2 ) 436 ENDIF 437 438 !---------------------------------------------------------------- 439 ! Vertical interpolation to the observation point 440 !---------------------------------------------------------------- 441 ista = prodatqc%npvsta(jobs,2) 442 iend = prodatqc%npvend(jobs,2) 443 CALL obs_int_z1d( kpk, & 444 & prodatqc%var(2)%mvk(ista:iend),& 445 & k1dint, iend - ista + 1, & 446 & prodatqc%var(2)%vdep(ista:iend),& 447 & zobsk, zobs2k, & 448 & prodatqc%var(2)%vmod(ista:iend),& 449 & pgdept, zobsmask2 ) 450 451 ENDIF 452 453 END DO 521 522 ! vertically interpolate all 4 corners 523 ista = prodatqc%npvsta(jobs,2) 524 iend = prodatqc%npvend(jobs,2) 525 inum_obs = iend - ista + 1 526 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 527 528 DO iin=1,2 529 DO ijn=1,2 530 531 IF ( k1dint == 1 ) THEN 532 CALL obs_int_z1d_spl( kpk, & 533 & zinm2(iin,ijn,:,iobs), & 534 & zobs2k, zgdept(iin,ijn,:,iobs), & 535 & zmask2(iin,ijn,:,iobs)) 536 ENDIF 537 538 CALL obs_level_search(kpk, & 539 & zgdept(iin,ijn,:,iobs), & 540 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 541 & iv_indic) 542 543 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 544 & prodatqc%var(2)%vdep(ista:iend), & 545 & zinm2(iin,ijn,:,iobs), & 546 & zobs2k, interp_corner(iin,ijn,:), & 547 & zgdept(iin,ijn,:,iobs), & 548 & zmask2(iin,ijn,:,iobs)) 549 550 ENDDO 551 ENDDO 552 553 ENDIF !idayend 554 555 ELSE 556 557 ! Point data 558 559 ! vertically interpolate all 4 corners 560 ista = prodatqc%npvsta(jobs,2) 561 iend = prodatqc%npvend(jobs,2) 562 inum_obs = iend - ista + 1 563 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 564 DO iin=1,2 565 DO ijn=1,2 566 567 IF ( k1dint == 1 ) THEN 568 CALL obs_int_z1d_spl( kpk, & 569 & zint2(iin,ijn,:,iobs),& 570 & zobs2k, zgdept(iin,ijn,:,iobs), & 571 & zmask2(iin,ijn,:,iobs)) 572 573 ENDIF 574 575 CALL obs_level_search(kpk, & 576 & zgdept(iin,ijn,:,iobs),& 577 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 578 & iv_indic) 579 580 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 581 & prodatqc%var(2)%vdep(ista:iend), & 582 & zint2(iin,ijn,:,iobs), & 583 & zobs2k,interp_corner(iin,ijn,:), & 584 & zgdept(iin,ijn,:,iobs), & 585 & zmask2(iin,ijn,:,iobs) ) 586 587 ENDDO 588 ENDDO 589 590 ENDIF 591 592 !------------------------------------------------------------- 593 ! Compute the horizontal interpolation for every profile level 594 !------------------------------------------------------------- 595 596 DO ikn=1,inum_obs 597 iend=ista+ikn-1 598 599 zweig(:,:,1) = 0._wp 600 601 ! This code forces the horizontal weights to be 602 ! zero IF the observation is below the bottom of the 603 ! corners of the interpolation nodes, Or if it is in 604 ! the mask. This is important for observations near 605 ! steep bathymetry 606 DO iin=1,2 607 DO ijn=1,2 608 609 depth_loop2: DO ik=kpk,2,-1 610 IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN 611 612 zweig(iin,ijn,1) = & 613 & zweig2(iin,ijn,1) * & 614 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 615 & - prodatqc%var(2)%vdep(iend)),0._wp) 616 617 EXIT depth_loop2 618 619 ENDIF 620 621 ENDDO depth_loop2 622 623 ENDDO 624 ENDDO 625 626 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 627 & prodatqc%var(2)%vmod(iend:iend) ) 628 629 ! Set QC flag for any observations found below the bottom 630 ! needed as the check here is more strict than that in obs_prep 631 IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 632 633 ENDDO 634 635 DEALLOCATE(interp_corner,iv_indic) 636 637 ENDIF 638 639 ENDDO 454 640 455 641 ! Deallocate the data for interpolation … … 466 652 & zmask2, & 467 653 & zint1, & 468 & zint2 & 654 & zint2, & 655 & zgdept, & 656 & zgdepw & 469 657 & ) 470 658 … … 481 669 END SUBROUTINE obs_prof_opt 482 670 483 SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 484 & ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, & 485 & kdailyavtypes ) 486 !!----------------------------------------------------------------------- 487 !! 488 !! *** ROUTINE obs_pro_opt *** 489 !! 490 !! ** Purpose : Compute the model counterpart of profiles 491 !! data by interpolating from the model grid to the 492 !! observation point. Generalised vertical coordinate version 493 !! 494 !! ** Method : Linearly interpolate to each observation point using 495 !! the model values at the corners of the surrounding grid box. 496 !! 497 !! First, model values on the model grid are interpolated vertically to the 498 !! Depths of the profile observations. Two vertical interpolation schemes are 499 !! available: 500 !! - linear (k1dint = 0) 501 !! - Cubic spline (k1dint = 1) 502 !! 503 !! 504 !! Secondly the interpolated values are interpolated horizontally to the 505 !! obs (lon, lat) point. 506 !! Several horizontal interpolation schemes are available: 507 !! - distance-weighted (great circle) (k2dint = 0) 508 !! - distance-weighted (small angle) (k2dint = 1) 509 !! - bilinear (geographical grid) (k2dint = 2) 510 !! - bilinear (quadrilateral grid) (k2dint = 3) 511 !! - polynomial (quadrilateral grid) (k2dint = 4) 512 !! 513 !! For the cubic spline the 2nd derivative of the interpolating 514 !! polynomial is computed before entering the vertical interpolation 515 !! routine. 516 !! 517 !! For ENACT moored buoy data (e.g., TAO), the model equivalent is 518 !! a daily mean model temperature field. So, we first compute 519 !! the mean, then interpolate only at the end of the day. 520 !! 521 !! This is the procedure to be used with generalised vertical model 522 !! coordinates (ie s-coordinates. It is ~4x slower than the equivalent 523 !! horizontal then vertical interpolation algorithm, but can deal with situations 524 !! where the model levels are not flat. 525 !! ONLY PERFORMED if ln_sco=.TRUE. 526 !! 527 !! Note: the in situ temperature observations must be converted 528 !! to potential temperature (the model variable) prior to 529 !! assimilation. 530 !!?????????????????????????????????????????????????????????????? 531 !! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR??? 532 !!?????????????????????????????????????????????????????????????? 533 !! 534 !! ** Action : 535 !! 536 !! History : 537 !! ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised 538 !! vertical coordinates 539 !!----------------------------------------------------------------------- 540 541 !! * Modules used 542 USE obs_profiles_def ! Definition of storage space for profile obs. 543 544 IMPLICIT NONE 545 546 !! * Arguments 547 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 548 INTEGER, INTENT(IN) :: kt ! Time step 549 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 550 INTEGER, INTENT(IN) :: kpj 551 INTEGER, INTENT(IN) :: kpk 552 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 553 ! (kit000-1 = restart time) 554 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 555 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 556 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 557 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 558 & ptn, & ! Model temperature field 559 & psn, & ! Model salinity field 560 & ptmask ! Land-sea mask 561 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 562 & pgdept, & ! Model array of depth T levels 563 & pgdepw ! Model array of depth W levels 564 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 565 & kdailyavtypes ! Types for daily averages 566 567 !! * Local declarations 568 INTEGER :: ji 569 INTEGER :: jj 570 INTEGER :: jk 571 INTEGER :: iico, ijco 572 INTEGER :: jobs 573 INTEGER :: inrc 574 INTEGER :: ipro 575 INTEGER :: idayend 576 INTEGER :: ista 577 INTEGER :: iend 578 INTEGER :: iobs 579 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 580 INTEGER, DIMENSION(imaxavtypes) :: & 581 & idailyavtypes 582 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 583 & igrdi, & 584 & igrdj 585 INTEGER :: & 586 & inum_obs 587 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 588 REAL(KIND=wp) :: zlam 589 REAL(KIND=wp) :: zphi 590 REAL(KIND=wp) :: zdaystp 591 REAL(KIND=wp), DIMENSION(kpk) :: & 592 & zobsmask, & 593 & zobsk, & 594 & zobs2k 595 REAL(KIND=wp), DIMENSION(2,2,1) :: & 596 & zweig, & 597 & l_zweig 598 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 599 & zmask, & 600 & zintt, & 601 & zints, & 602 & zinmt, & 603 & zgdept,& 604 & zgdepw,& 605 & zinms 606 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 607 & zglam, & 608 & zgphi 609 REAL(KIND=wp), DIMENSION(1) :: zmsk_1 610 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 611 612 !------------------------------------------------------------------------ 613 ! Local initialization 614 !------------------------------------------------------------------------ 615 ! ... Record and data counters 616 inrc = kt - kit000 + 2 617 ipro = prodatqc%npstp(inrc) 618 619 ! Daily average types 620 IF ( PRESENT(kdailyavtypes) ) THEN 621 idailyavtypes(:) = kdailyavtypes(:) 622 ELSE 623 idailyavtypes(:) = -1 624 ENDIF 625 626 ! Initialize daily mean for first time-step 627 idayend = MOD( kt - kit000 + 1, kdaystp ) 628 629 ! Added kt == 0 test to catch restart case 630 IF ( idayend == 1 .OR. kt == 0) THEN 631 632 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 633 DO jk = 1, jpk 634 DO jj = 1, jpj 635 DO ji = 1, jpi 636 prodatqc%vdmean(ji,jj,jk,1) = 0.0 637 prodatqc%vdmean(ji,jj,jk,2) = 0.0 638 END DO 639 END DO 640 END DO 641 642 ENDIF 643 644 DO jk = 1, jpk 645 DO jj = 1, jpj 646 DO ji = 1, jpi 647 ! Increment the temperature field for computing daily mean 648 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 649 & + ptn(ji,jj,jk) 650 ! Increment the salinity field for computing daily mean 651 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 652 & + psn(ji,jj,jk) 653 END DO 654 END DO 655 END DO 656 657 ! Compute the daily mean at the end of day 658 zdaystp = 1.0 / REAL( kdaystp ) 659 IF ( idayend == 0 ) THEN 660 DO jk = 1, jpk 661 DO jj = 1, jpj 662 DO ji = 1, jpi 663 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 664 & * zdaystp 665 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 666 & * zdaystp 667 END DO 668 END DO 669 END DO 670 ENDIF 671 672 ! Get the data for interpolation 673 ALLOCATE( & 674 & igrdi(2,2,ipro), & 675 & igrdj(2,2,ipro), & 676 & zglam(2,2,ipro), & 677 & zgphi(2,2,ipro), & 678 & zmask(2,2,kpk,ipro), & 679 & zintt(2,2,kpk,ipro), & 680 & zints(2,2,kpk,ipro), & 681 & zgdept(2,2,kpk,ipro), & 682 & zgdepw(2,2,kpk,ipro) & 683 & ) 684 685 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 686 iobs = jobs - prodatqc%nprofup 687 igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 688 igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 689 igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 690 igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 691 igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 692 igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 693 igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 694 igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 695 END DO 696 697 ! Initialise depth arrays 698 zgdept = 0.0 699 zgdepw = 0.0 700 701 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, glamt, zglam ) 702 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, gphit, zgphi ) 703 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptmask,zmask ) 704 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptn, zintt ) 705 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, psn, zints ) 706 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept(:,:,:), & 707 & zgdept ) 708 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw(:,:,:), & 709 & zgdepw ) 710 711 ! At the end of the day also get interpolated means 712 IF ( idayend == 0 ) THEN 713 714 ALLOCATE( & 715 & zinmt(2,2,kpk,ipro), & 716 & zinms(2,2,kpk,ipro) & 717 & ) 718 719 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 720 & prodatqc%vdmean(:,:,:,1), zinmt ) 721 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 722 & prodatqc%vdmean(:,:,:,2), zinms ) 723 724 ENDIF 725 726 ! Return if no observations to process 727 ! Has to be done after comm commands to ensure processors 728 ! stay in sync 729 IF ( ipro == 0 ) RETURN 730 731 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 732 733 iobs = jobs - prodatqc%nprofup 734 735 IF ( kt /= prodatqc%mstp(jobs) ) THEN 736 737 IF(lwp) THEN 738 WRITE(numout,*) 739 WRITE(numout,*) ' E R R O R : Observation', & 740 & ' time step is not consistent with the', & 741 & ' model time step' 742 WRITE(numout,*) ' =========' 743 WRITE(numout,*) 744 WRITE(numout,*) ' Record = ', jobs, & 745 & ' kt = ', kt, & 746 & ' mstp = ', prodatqc%mstp(jobs), & 747 & ' ntyp = ', prodatqc%ntyp(jobs) 748 ENDIF 749 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 750 ENDIF 751 752 zlam = prodatqc%rlam(jobs) 753 zphi = prodatqc%rphi(jobs) 754 755 ! Horizontal weights 756 ! Only calculated once, for both T and S. 757 ! Masked values are calculated later. 758 759 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 760 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 761 762 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 763 & zglam(:,:,iobs), zgphi(:,:,iobs), & 764 & zmask(:,:,1,iobs), zweig, zmsk_1 ) 765 766 ENDIF 767 768 ! IF zmsk_1 = 0; then ob is on land 769 IF (zmsk_1(1) < 0.1) THEN 770 WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask' 771 772 ELSE 773 774 ! Temperature 775 776 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 777 778 zobsk(:) = obfillflt 779 780 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 781 782 IF ( idayend == 0 ) THEN 783 784 ! Daily averaged moored buoy (MRB) data 785 786 ! vertically interpolate all 4 corners 787 ista = prodatqc%npvsta(jobs,1) 788 iend = prodatqc%npvend(jobs,1) 789 inum_obs = iend - ista + 1 790 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 791 792 DO iin=1,2 793 DO ijn=1,2 794 795 796 797 IF ( k1dint == 1 ) THEN 798 CALL obs_int_z1d_spl( kpk, & 799 & zinmt(iin,ijn,:,iobs), & 800 & zobs2k, zgdept(iin,ijn,:,iobs), & 801 & zmask(iin,ijn,:,iobs)) 802 ENDIF 803 804 CALL obs_level_search(kpk, & 805 & zgdept(iin,ijn,:,iobs), & 806 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 807 & iv_indic) 808 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 809 & prodatqc%var(1)%vdep(ista:iend), & 810 & zinmt(iin,ijn,:,iobs), & 811 & zobs2k, interp_corner(iin,ijn,:), & 812 & zgdept(iin,ijn,:,iobs), & 813 & zmask(iin,ijn,:,iobs)) 814 815 ENDDO 816 ENDDO 817 818 819 ELSE 820 821 CALL ctl_stop( ' A nonzero' // & 822 & ' number of profile T BUOY data should' // & 823 & ' only occur at the end of a given day' ) 824 825 ENDIF 826 827 ELSE 828 829 ! Point data 830 831 ! vertically interpolate all 4 corners 832 ista = prodatqc%npvsta(jobs,1) 833 iend = prodatqc%npvend(jobs,1) 834 inum_obs = iend - ista + 1 835 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 836 DO iin=1,2 837 DO ijn=1,2 838 839 840 IF ( k1dint == 1 ) THEN 841 CALL obs_int_z1d_spl( kpk, & 842 & zintt(iin,ijn,:,iobs),& 843 & zobs2k, zgdept(iin,ijn,:,iobs), & 844 & zmask(iin,ijn,:,iobs)) 845 846 ENDIF 847 848 CALL obs_level_search(kpk, & 849 & zgdept(iin,ijn,:,iobs),& 850 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 851 & iv_indic) 852 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 853 & prodatqc%var(1)%vdep(ista:iend), & 854 & zintt(iin,ijn,:,iobs), & 855 & zobs2k,interp_corner(iin,ijn,:), & 856 & zgdept(iin,ijn,:,iobs), & 857 & zmask(iin,ijn,:,iobs) ) 858 859 ENDDO 860 ENDDO 861 862 ENDIF 863 864 !------------------------------------------------------------- 865 ! Compute the horizontal interpolation for every profile level 866 !------------------------------------------------------------- 867 868 DO ikn=1,inum_obs 869 iend=ista+ikn-1 870 871 l_zweig(:,:,1) = 0._wp 872 873 ! This code forces the horizontal weights to be 874 ! zero IF the observation is below the bottom of the 875 ! corners of the interpolation nodes, Or if it is in 876 ! the mask. This is important for observations are near 877 ! steep bathymetry 878 DO iin=1,2 879 DO ijn=1,2 880 881 depth_loop1: DO ik=kpk,2,-1 882 IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN 883 884 l_zweig(iin,ijn,1) = & 885 & zweig(iin,ijn,1) * & 886 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 887 & - prodatqc%var(1)%vdep(iend)),0._wp) 888 889 EXIT depth_loop1 890 ENDIF 891 ENDDO depth_loop1 892 893 ENDDO 894 ENDDO 895 896 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 897 & prodatqc%var(1)%vmod(iend:iend) ) 898 899 ENDDO 900 901 902 DEALLOCATE(interp_corner,iv_indic) 903 904 ENDIF 905 906 907 ! Salinity 908 909 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 910 911 zobsk(:) = obfillflt 912 913 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 914 915 IF ( idayend == 0 ) THEN 916 917 ! Daily averaged moored buoy (MRB) data 918 919 ! vertically interpolate all 4 corners 920 ista = prodatqc%npvsta(iobs,2) 921 iend = prodatqc%npvend(iobs,2) 922 inum_obs = iend - ista + 1 923 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 924 925 DO iin=1,2 926 DO ijn=1,2 927 928 929 930 IF ( k1dint == 1 ) THEN 931 CALL obs_int_z1d_spl( kpk, & 932 & zinms(iin,ijn,:,iobs), & 933 & zobs2k, zgdept(iin,ijn,:,iobs), & 934 & zmask(iin,ijn,:,iobs)) 935 ENDIF 936 937 CALL obs_level_search(kpk, & 938 & zgdept(iin,ijn,:,iobs), & 939 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 940 & iv_indic) 941 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 942 & prodatqc%var(2)%vdep(ista:iend), & 943 & zinms(iin,ijn,:,iobs), & 944 & zobs2k, interp_corner(iin,ijn,:), & 945 & zgdept(iin,ijn,:,iobs), & 946 & zmask(iin,ijn,:,iobs)) 947 948 ENDDO 949 ENDDO 950 951 952 ELSE 953 954 CALL ctl_stop( ' A nonzero' // & 955 & ' number of profile T BUOY data should' // & 956 & ' only occur at the end of a given day' ) 957 958 ENDIF 959 960 ELSE 961 962 ! Point data 963 964 ! vertically interpolate all 4 corners 965 ista = prodatqc%npvsta(jobs,2) 966 iend = prodatqc%npvend(jobs,2) 967 inum_obs = iend - ista + 1 968 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 969 970 DO iin=1,2 971 DO ijn=1,2 972 973 974 IF ( k1dint == 1 ) THEN 975 CALL obs_int_z1d_spl( kpk, & 976 & zints(iin,ijn,:,iobs),& 977 & zobs2k, zgdept(iin,ijn,:,iobs), & 978 & zmask(iin,ijn,:,iobs)) 979 980 ENDIF 981 982 CALL obs_level_search(kpk, & 983 & zgdept(iin,ijn,:,iobs),& 984 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 985 & iv_indic) 986 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 987 & prodatqc%var(2)%vdep(ista:iend), & 988 & zints(iin,ijn,:,iobs), & 989 & zobs2k,interp_corner(iin,ijn,:), & 990 & zgdept(iin,ijn,:,iobs), & 991 & zmask(iin,ijn,:,iobs) ) 992 993 ENDDO 994 ENDDO 995 996 ENDIF 997 998 !------------------------------------------------------------- 999 ! Compute the horizontal interpolation for every profile level 1000 !------------------------------------------------------------- 1001 1002 DO ikn=1,inum_obs 1003 iend=ista+ikn-1 1004 1005 l_zweig(:,:,1) = 0._wp 1006 1007 ! This code forces the horizontal weights to be 1008 ! zero IF the observation is below the bottom of the 1009 ! corners of the interpolation nodes, Or if it is in 1010 ! the mask. This is important for observations are near 1011 ! steep bathymetry 1012 DO iin=1,2 1013 DO ijn=1,2 1014 1015 depth_loop2: DO ik=kpk,2,-1 1016 IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN 1017 1018 l_zweig(iin,ijn,1) = & 1019 & zweig(iin,ijn,1) * & 1020 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 1021 & - prodatqc%var(2)%vdep(iend)),0._wp) 1022 1023 EXIT depth_loop2 1024 ENDIF 1025 ENDDO depth_loop2 1026 1027 ENDDO 1028 ENDDO 1029 1030 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 1031 & prodatqc%var(2)%vmod(iend:iend) ) 1032 1033 ENDDO 1034 1035 1036 DEALLOCATE(interp_corner,iv_indic) 1037 1038 ENDIF 1039 1040 ENDIF 1041 1042 END DO 1043 1044 ! Deallocate the data for interpolation 1045 DEALLOCATE( & 1046 & igrdi, & 1047 & igrdj, & 1048 & zglam, & 1049 & zgphi, & 1050 & zmask, & 1051 & zintt, & 1052 & zints, & 1053 & zgdept,& 1054 & zgdepw & 1055 & ) 1056 ! At the end of the day also get interpolated means 1057 IF ( idayend == 0 ) THEN 1058 DEALLOCATE( & 1059 & zinmt, & 1060 & zinms & 1061 & ) 1062 ENDIF 1063 1064 prodatqc%nprofup = prodatqc%nprofup + ipro 1065 1066 END SUBROUTINE obs_pro_sco_opt 1067 1068 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 1069 & kit000, kdaystp, psurf, psurfmask, & 1070 & k2dint, ldnightav ) 671 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 672 & kit000, kdaystp, psurf, psurfmask, & 673 & k2dint, ldnightav, plamscl, pphiscl, & 674 & lindegrees ) 1071 675 1072 676 !!----------------------------------------------------------------------- … … 1090 694 !! - polynomial (quadrilateral grid) (k2dint = 4) 1091 695 !! 696 !! Two horizontal averaging schemes are also available: 697 !! - weighted radial footprint (k2dint = 5) 698 !! - weighted rectangular footprint (k2dint = 6) 699 !! 1092 700 !! 1093 701 !! ** Action : … … 1096 704 !! ! 07-03 (A. Weaver) 1097 705 !! ! 15-02 (M. Martin) Combined routine for surface types 706 !! ! 17-03 (M. Martin) Added horizontal averaging options 1098 707 !!----------------------------------------------------------------------- 1099 708 … … 1117 726 & psurfmask ! Land-sea mask 1118 727 LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 728 REAL(KIND=wp), INTENT(IN) :: & 729 & plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions 730 & pphiscl ! This is the full width (rather than half-width) 731 LOGICAL, INTENT(IN) :: & 732 & lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 1119 733 1120 734 !! * Local declarations … … 1125 739 INTEGER :: isurf 1126 740 INTEGER :: iobs 741 INTEGER :: imaxifp, imaxjfp 742 INTEGER :: imodi, imodj 1127 743 INTEGER :: idayend 1128 744 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 1129 & igrdi, & 1130 & igrdj 745 & igrdi, & 746 & igrdj, & 747 & igrdip1, & 748 & igrdjp1 1131 749 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1132 750 & icount_night, & … … 1136 754 REAL(wp), DIMENSION(1) :: zext, zobsmask 1137 755 REAL(wp) :: zdaystp 1138 REAL(wp), DIMENSION(2,2,1) :: &1139 & zweig1140 756 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 757 & zweig, & 1141 758 & zmask, & 1142 759 & zsurf, & 1143 760 & zsurfm, & 761 & zsurftmp, & 1144 762 & zglam, & 1145 & zgphi 763 & zgphi, & 764 & zglamf, & 765 & zgphif 766 1146 767 REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1147 768 & zintmp, & … … 1155 776 inrc = kt - kit000 + 2 1156 777 isurf = surfdataqc%nsstp(inrc) 778 779 ! Work out the maximum footprint size for the 780 ! interpolation/averaging in model grid-points - has to be even. 781 782 CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 783 1157 784 1158 785 IF ( ldnightav ) THEN … … 1221 848 1222 849 ALLOCATE( & 1223 & igrdi(2,2,isurf), & 1224 & igrdj(2,2,isurf), & 1225 & zglam(2,2,isurf), & 1226 & zgphi(2,2,isurf), & 1227 & zmask(2,2,isurf), & 1228 & zsurf(2,2,isurf) & 850 & zweig(imaxifp,imaxjfp,1), & 851 & igrdi(imaxifp,imaxjfp,isurf), & 852 & igrdj(imaxifp,imaxjfp,isurf), & 853 & zglam(imaxifp,imaxjfp,isurf), & 854 & zgphi(imaxifp,imaxjfp,isurf), & 855 & zmask(imaxifp,imaxjfp,isurf), & 856 & zsurf(imaxifp,imaxjfp,isurf), & 857 & zsurftmp(imaxifp,imaxjfp,isurf), & 858 & zglamf(imaxifp+1,imaxjfp+1,isurf), & 859 & zgphif(imaxifp+1,imaxjfp+1,isurf), & 860 & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 861 & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 1229 862 & ) 1230 863 1231 864 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 1232 865 iobs = jobs - surfdataqc%nsurfup 1233 igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1 1234 igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1 1235 igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1 1236 igrdj(1,2,iobs) = surfdataqc%mj(jobs) 1237 igrdi(2,1,iobs) = surfdataqc%mi(jobs) 1238 igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1 1239 igrdi(2,2,iobs) = surfdataqc%mi(jobs) 1240 igrdj(2,2,iobs) = surfdataqc%mj(jobs) 866 DO ji = 0, imaxifp 867 imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 868 869 !Deal with wrap around in longitude 870 IF ( imodi < 1 ) imodi = imodi + jpiglo 871 IF ( imodi > jpiglo ) imodi = imodi - jpiglo 872 873 DO jj = 0, imaxjfp 874 imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 875 !If model values are out of the domain to the north/south then 876 !set them to be the edge of the domain 877 IF ( imodj < 1 ) imodj = 1 878 IF ( imodj > jpjglo ) imodj = jpjglo 879 880 igrdip1(ji+1,jj+1,iobs) = imodi 881 igrdjp1(ji+1,jj+1,iobs) = imodj 882 883 IF ( ji >= 1 .AND. jj >= 1 ) THEN 884 igrdi(ji,jj,iobs) = imodi 885 igrdj(ji,jj,iobs) = imodj 886 ENDIF 887 888 END DO 889 END DO 1241 890 END DO 1242 891 1243 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &892 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1244 893 & igrdi, igrdj, glamt, zglam ) 1245 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &894 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1246 895 & igrdi, igrdj, gphit, zgphi ) 1247 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &896 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1248 897 & igrdi, igrdj, psurfmask, zmask ) 1249 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &898 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1250 899 & igrdi, igrdj, psurf, zsurf ) 900 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 901 & igrdip1, igrdjp1, glamf, zglamf ) 902 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 903 & igrdip1, igrdjp1, gphif, zgphif ) 1251 904 1252 905 ! At the end of the day get interpolated means 1253 IF (ldnightav ) THEN 1254 IF ( idayend == 0 ) THEN 1255 1256 ALLOCATE( & 1257 & zsurfm(2,2,isurf) & 1258 & ) 1259 1260 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, & 1261 & surfdataqc%vdmean(:,:), zsurfm ) 1262 1263 ENDIF 906 IF ( idayend == 0 .AND. ldnightav ) THEN 907 908 ALLOCATE( & 909 & zsurfm(imaxifp,imaxjfp,isurf) & 910 & ) 911 912 CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 913 & surfdataqc%vdmean(:,:), zsurfm ) 914 1264 915 ENDIF 1265 916 … … 1290 941 zphi = surfdataqc%rphi(jobs) 1291 942 1292 ! Get weights to interpolate the model value to the observation point1293 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &1294 & zglam(:,:,iobs), zgphi(:,:,iobs), &1295 & zmask(:,:,iobs), zweig, zobsmask )1296 1297 ! Interpolate the model field to the observation point1298 943 IF ( ldnightav .AND. idayend == 0 ) THEN 1299 944 ! Night-time averaged data 1300 CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext)945 zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 1301 946 ELSE 1302 CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs), zext ) 947 zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 948 ENDIF 949 950 IF ( k2dint <= 4 ) THEN 951 952 ! Get weights to interpolate the model value to the observation point 953 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 954 & zglam(:,:,iobs), zgphi(:,:,iobs), & 955 & zmask(:,:,iobs), zweig, zobsmask ) 956 957 ! Interpolate the model value to the observation point 958 CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 959 960 ELSE 961 962 ! Get weights to average the model SLA to the observation footprint 963 CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, & 964 & zglam(:,:,iobs), zgphi(:,:,iobs), & 965 & zglamf(:,:,iobs), zgphif(:,:,iobs), & 966 & zmask(:,:,iobs), plamscl, pphiscl, & 967 & lindegrees, zweig, zobsmask ) 968 969 ! Average the model SST to the observation footprint 970 CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 971 & zweig, zsurftmp(:,:,iobs), zext ) 972 1303 973 ENDIF 1304 974 … … 1310 980 surfdataqc%rmod(jobs,1) = zext(1) 1311 981 ENDIF 982 983 IF ( zext(1) == obfillflt ) THEN 984 ! If the observation value is a fill value, set QC flag to bad 985 surfdataqc%nqc(jobs) = 4 986 ENDIF 1312 987 1313 988 END DO … … 1315 990 ! Deallocate the data for interpolation 1316 991 DEALLOCATE( & 992 & zweig, & 1317 993 & igrdi, & 1318 994 & igrdj, & … … 1320 996 & zgphi, & 1321 997 & zmask, & 1322 & zsurf & 998 & zsurf, & 999 & zsurftmp, & 1000 & zglamf, & 1001 & zgphif, & 1002 & igrdip1,& 1003 & igrdjp1 & 1323 1004 & ) 1324 1005 1325 1006 ! At the end of the day also deallocate night-time mean array 1326 IF ( ldnightav ) THEN 1327 IF ( idayend == 0 ) THEN 1328 DEALLOCATE( & 1329 & zsurfm & 1330 & ) 1331 ENDIF 1007 IF ( idayend == 0 .AND. ldnightav ) THEN 1008 DEALLOCATE( & 1009 & zsurfm & 1010 & ) 1332 1011 ENDIF 1333 1012 -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
r7646 r8992 23 23 USE obs_oper ! Observation operators 24 24 USE lib_mpp, ONLY : ctl_warn, ctl_stop 25 USE bdy_oce, ONLY : & ! Boundary information 26 idx_bdy, nb_bdy, ln_bdy 25 27 26 28 IMPLICIT NONE … … 40 42 CONTAINS 41 43 42 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 44 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 45 kqc_cutoff ) 43 46 !!---------------------------------------------------------------------- 44 47 !! *** ROUTINE obs_pre_sla *** … … 57 60 !! ! 2015-02 (M. Martin) Combined routine for surface types. 58 61 !!---------------------------------------------------------------------- 62 !! * Modules used 59 63 USE par_oce ! Ocean parameters 60 64 USE dom_oce, ONLY : glamt, gphit, tmask, nproc ! Geographical information … … 63 67 TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of surface data not failing screening 64 68 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 65 ! 69 LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting obs near the boundary 70 INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value 71 !! * Local declarations 72 INTEGER :: iqc_cutoff = 255 ! cut off for QC value 66 73 INTEGER :: iyea0 ! Initial date 67 74 INTEGER :: imon0 ! - (year, month, day, hour, minute) … … 76 83 INTEGER :: inlasobs ! - close to land 77 84 INTEGER :: igrdobs ! - fail the grid search 85 INTEGER :: ibdysobs ! - close to open boundary 78 86 ! Global counters for observations that 79 87 INTEGER :: iotdobsmpp ! - outside time domain … … 82 90 INTEGER :: inlasobsmpp ! - close to land 83 91 INTEGER :: igrdobsmpp ! - fail the grid search 84 LOGICAL, DIMENSION(:), ALLOCATABLE :: llvalid ! SLA data selection 92 INTEGER :: ibdysobsmpp ! - close to open boundary 93 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 94 & llvalid ! SLA data selection 85 95 INTEGER :: jobs ! Obs. loop variable 86 96 INTEGER :: jstp ! Time loop variable … … 107 117 ilansobs = 0 108 118 inlasobs = 0 119 ibdysobs = 0 120 121 ! Set QC cutoff to optional value if provided 122 IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 109 123 110 124 ! ----------------------------------------------------------------------- … … 140 154 & tmask(:,:,1), surfdata%nqc, & 141 155 & iosdsobs, ilansobs, & 142 & inlasobs, ld_nea ) 156 & inlasobs, ld_nea, & 157 & ibdysobs, ld_bound_reject, & 158 & iqc_cutoff ) 143 159 144 160 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 145 161 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 146 162 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 163 CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 147 164 148 165 ! ----------------------------------------------------------------------- … … 155 172 ALLOCATE( llvalid(surfdata%nsurf) ) 156 173 157 ! We want all data which has qc flags <= 10158 159 llvalid(:) = ( surfdata%nqc(:) <= 10)174 ! We want all data which has qc flags <= iqc_cutoff 175 176 llvalid(:) = ( surfdata%nqc(:) <= iqc_cutoff ) 160 177 161 178 ! The actual copying … … 190 207 & inlasobsmpp 191 208 ENDIF 209 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 210 & ibdysobsmpp 192 211 WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted = ', & 193 212 & surfdataqc%nsurfmpp … … 225 244 & kpi, kpj, kpk, & 226 245 & zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2, & 227 & ld_nea, kdailyavtypes)246 & ld_nea, ld_bound_reject, kdailyavtypes, kqc_cutoff ) 228 247 229 248 !!---------------------------------------------------------------------- … … 241 260 !! 242 261 !!---------------------------------------------------------------------- 243 USE par_oce ! Ocean parameters 244 USE dom_oce, ONLY : gdept_1d, nproc ! Geographical information 262 !! * Modules used 263 USE par_oce ! Ocean parameters 264 USE dom_oce, ONLY : & ! Geographical information 265 & gdept_1d, & 266 & nproc 245 267 246 268 !! * Arguments … … 250 272 LOGICAL, INTENT(IN) :: ld_var2 251 273 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 274 LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting observations near the boundary 252 275 INTEGER, INTENT(IN) :: kpi, kpj, kpk ! Local domain sizes 253 276 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & … … 261 284 & pgphi1, & 262 285 & pgphi2 286 INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value 263 287 264 288 !! * Local declarations 289 INTEGER :: iqc_cutoff = 255 ! cut off for QC value 265 290 INTEGER :: iyea0 ! Initial date 266 291 INTEGER :: imon0 ! - (year, month, day, hour, minute) … … 277 302 INTEGER :: inlav1obs ! - close to land (variable 1) 278 303 INTEGER :: inlav2obs ! - close to land (variable 2) 304 INTEGER :: ibdyv1obs ! - boundary (variable 1) 305 INTEGER :: ibdyv2obs ! - boundary (variable 2) 279 306 INTEGER :: igrdobs ! - fail the grid search 280 307 INTEGER :: iuvchku ! - reject u if v rejected and vice versa … … 288 315 INTEGER :: inlav1obsmpp ! - close to land (variable 1) 289 316 INTEGER :: inlav2obsmpp ! - close to land (variable 2) 317 INTEGER :: ibdyv1obsmpp ! - boundary (variable 1) 318 INTEGER :: ibdyv2obsmpp ! - boundary (variable 2) 290 319 INTEGER :: igrdobsmpp ! - fail the grid search 291 320 INTEGER :: iuvchkumpp ! - reject var1 if var2 rejected and vice versa … … 322 351 inlav1obs = 0 323 352 inlav2obs = 0 324 iuvchku = 0 325 iuvchkv = 0 353 ibdyv1obs = 0 354 ibdyv2obs = 0 355 iuvchku = 0 356 iuvchkv = 0 357 358 359 ! Set QC cutoff to optional value if provided 360 IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 326 361 327 362 ! ----------------------------------------------------------------------- … … 335 370 & profdata%nday, profdata%nhou, profdata%nmin, & 336 371 & profdata%ntyp, profdata%nqc, profdata%mstp, & 337 & iotdobs, kdailyavtypes = kdailyavtypes ) 372 & iotdobs, kdailyavtypes = kdailyavtypes, & 373 & kqc_cutoff = iqc_cutoff ) 338 374 ELSE 339 375 CALL obs_coo_tim_prof( icycle, & … … 342 378 & profdata%nday, profdata%nhou, profdata%nmin, & 343 379 & profdata%ntyp, profdata%nqc, profdata%mstp, & 344 & iotdobs )380 & iotdobs, kqc_cutoff = iqc_cutoff ) 345 381 ENDIF 346 382 … … 359 395 360 396 ! ----------------------------------------------------------------------- 361 ! Reject all observations for profiles with nqc > 10362 ! ----------------------------------------------------------------------- 363 364 CALL obs_pro_rej( profdata )397 ! Reject all observations for profiles with nqc > iqc_cutoff 398 ! ----------------------------------------------------------------------- 399 400 CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff ) 365 401 366 402 ! ----------------------------------------------------------------------- … … 381 417 & gdept_1d, zmask1, & 382 418 & profdata%nqc, profdata%var(1)%nvqc, & 383 & iosdv1obs, ilanv1obs, & 384 & inlav1obs, ld_nea ) 419 & iosdv1obs, ilanv1obs, & 420 & inlav1obs, ld_nea, & 421 & ibdyv1obs, ld_bound_reject, & 422 & iqc_cutoff ) 385 423 386 424 CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 387 425 CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 388 426 CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 427 CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp ) 389 428 390 429 ! Variable 2 … … 400 439 & gdept_1d, zmask2, & 401 440 & profdata%nqc, profdata%var(2)%nvqc, & 402 & iosdv2obs, ilanv2obs, & 403 & inlav2obs, ld_nea ) 441 & iosdv2obs, ilanv2obs, & 442 & inlav2obs, ld_nea, & 443 & ibdyv2obs, ld_bound_reject, & 444 & iqc_cutoff ) 404 445 405 446 CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 406 447 CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 407 448 CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 449 CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 408 450 409 451 ! ----------------------------------------------------------------------- … … 412 454 413 455 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 414 CALL obs_uv_rej( profdata, iuvchku, iuvchkv )456 CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 415 457 CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 416 458 CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) … … 429 471 END DO 430 472 431 ! We want all data which has qc flags = 0432 433 llvalid%luse(:) = ( profdata%nqc(:) <= 10)473 ! We want all data which has qc flags <= iqc_cutoff 474 475 llvalid%luse(:) = ( profdata%nqc(:) <= iqc_cutoff ) 434 476 DO jvar = 1,profdata%nvar 435 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10)477 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff ) 436 478 END DO 437 479 … … 475 517 & iuvchku 476 518 ENDIF 519 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',& 520 & ibdyv1obsmpp 477 521 WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted = ', & 478 522 & prodatqc%nvprotmpp(1) … … 492 536 & iuvchkv 493 537 ENDIF 538 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 539 & ibdyv2obsmpp 494 540 WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted = ', & 495 541 & prodatqc%nvprotmpp(2) … … 644 690 & .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN 645 691 kobsstp(jobs) = -1 646 kobsqc(jobs) = kobsqc(jobs) + 11692 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 647 693 kotdobs = kotdobs + 1 648 694 CYCLE … … 695 741 IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) & 696 742 & .OR.( kobsstp(jobs) > nitend ) ) THEN 697 kobsqc(jobs) = kobsqc(jobs) + 12743 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 698 744 kotdobs = kotdobs + 1 699 745 CYCLE … … 739 785 & kobsno, & 740 786 & kobsyea, kobsmon, kobsday, kobshou, kobsmin, & 741 & ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes ) 787 & ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, & 788 & kqc_cutoff ) 742 789 !!---------------------------------------------------------------------- 743 790 !! *** ROUTINE obs_coo_tim *** … … 783 830 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 784 831 & kdailyavtypes ! Types for daily averages 832 INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff ! QC cutoff value 833 785 834 !! * Local declarations 786 835 INTEGER :: jobs 836 INTEGER :: iqc_cutoff=255 787 837 788 838 !----------------------------------------------------------------------- … … 803 853 DO jobs = 1, kobsno 804 854 805 IF ( kobsqc(jobs) <= 10) THEN855 IF ( kobsqc(jobs) <= iqc_cutoff ) THEN 806 856 807 857 IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.& 808 858 & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN 809 kobsqc(jobs) = kobsqc(jobs) + 14859 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 810 860 kotdobs = kotdobs + 1 811 861 CYCLE … … 850 900 DO jobs = 1, kobsno 851 901 IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN 852 kobsqc(jobs) = kobsqc(jobs) + 18902 kobsqc(jobs) = IBSET(kobsqc(jobs),12) 853 903 kgrdobs = kgrdobs + 1 854 904 ENDIF … … 861 911 & plam, pphi, pmask, & 862 912 & kobsqc, kosdobs, klanobs, & 863 & knlaobs,ld_nea ) 913 & knlaobs,ld_nea, & 914 & kbdyobs,ld_bound_reject, & 915 & kqc_cutoff ) 864 916 !!---------------------------------------------------------------------- 865 917 !! *** ROUTINE obs_coo_spc_2d *** … … 894 946 INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 895 947 & kobsqc ! Observation quality control 896 INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain 897 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 898 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 899 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 948 INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain 949 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 950 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 951 INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary 952 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 953 LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary 954 INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value 955 900 956 !! * Local declarations 901 957 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 902 958 & zgmsk ! Grid mask 959 960 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 961 & zbmsk ! Boundary mask 962 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 903 963 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 904 964 & zglam, & ! Model longitude at grid points … … 917 977 ! For invalid points use 2,2 918 978 919 IF ( kobsqc(jobs) >= 10) THEN979 IF ( kobsqc(jobs) >= kqc_cutoff ) THEN 920 980 921 981 igrdi(1,1,jobs) = 1 … … 942 1002 943 1003 END DO 1004 1005 IF (ln_bdy) THEN 1006 ! Create a mask grid points in boundary rim 1007 IF (ld_bound_reject) THEN 1008 zbdymask(:,:) = 1.0_wp 1009 DO ji = 1, nb_bdy 1010 DO jj = 1, idx_bdy(ji)%nblen(1) 1011 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 1012 ENDDO 1013 ENDDO 1014 1015 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 1016 ENDIF 1017 ENDIF 1018 944 1019 945 1020 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) … … 950 1025 951 1026 ! Skip bad observations 952 IF ( kobsqc(jobs) >= 10) CYCLE1027 IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE 953 1028 954 1029 ! Flag if the observation falls outside the model spatial domain … … 957 1032 & .OR. ( pobsphi(jobs) < -90. ) & 958 1033 & .OR. ( pobsphi(jobs) > 90. ) ) THEN 959 kobsqc(jobs) = kobsqc(jobs) + 111034 kobsqc(jobs) = IBSET(kobsqc(jobs),11) 960 1035 kosdobs = kosdobs + 1 961 1036 CYCLE … … 964 1039 ! Flag if the observation falls with a model land cell 965 1040 IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 966 kobsqc(jobs) = kobsqc(jobs) + 121041 kobsqc(jobs) = IBSET(kobsqc(jobs),10) 967 1042 klanobs = klanobs + 1 968 1043 CYCLE … … 978 1053 IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 979 1054 & .AND. & 980 & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp )&981 & ) THEN1055 & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) & 1056 & < 1.0e-6_wp ) ) THEN 982 1057 lgridobs = .TRUE. 983 1058 iig = ji … … 992 1067 IF (lgridobs) THEN 993 1068 IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 994 kobsqc(jobs) = kobsqc(jobs) + 121069 kobsqc(jobs) = IBSET(kobsqc(jobs),10) 995 1070 klanobs = klanobs + 1 996 1071 CYCLE … … 1000 1075 ! Flag if the observation falls is close to land 1001 1076 IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN 1002 IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 141003 1077 knlaobs = knlaobs + 1 1004 CYCLE 1078 IF (ld_nea) THEN 1079 kobsqc(jobs) = IBSET(kobsqc(jobs),9) 1080 CYCLE 1081 ENDIF 1082 ENDIF 1083 1084 IF (ln_bdy) THEN 1085 ! Flag if the observation falls close to the boundary rim 1086 IF (ld_bound_reject) THEN 1087 IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 1088 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 1089 kbdyobs = kbdyobs + 1 1090 CYCLE 1091 ENDIF 1092 ! for observations on the grid... 1093 IF (lgridobs) THEN 1094 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 1095 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 1096 kbdyobs = kbdyobs + 1 1097 CYCLE 1098 ENDIF 1099 ENDIF 1100 ENDIF 1005 1101 ENDIF 1006 1102 … … 1015 1111 & plam, pphi, pdep, pmask, & 1016 1112 & kpobsqc, kobsqc, kosdobs, & 1017 & klanobs, knlaobs, ld_nea ) 1113 & klanobs, knlaobs, ld_nea, & 1114 & kbdyobs, ld_bound_reject, & 1115 & kqc_cutoff ) 1018 1116 !!---------------------------------------------------------------------- 1019 1117 !! *** ROUTINE obs_coo_spc_3d *** … … 1077 1175 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 1078 1176 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 1177 INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary 1079 1178 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 1179 LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary 1180 INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value 1181 1080 1182 !! * Local declarations 1081 1183 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1082 1184 & zgmsk ! Grid mask 1185 REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 1186 & zbmsk ! Boundary mask 1187 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 1083 1188 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1084 1189 & zgdepw … … 1100 1205 ! For invalid points use 2,2 1101 1206 1102 IF ( kpobsqc(jobs) >= 10) THEN1207 IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN 1103 1208 1104 1209 igrdi(1,1,jobs) = 1 … … 1125 1230 1126 1231 END DO 1232 1233 IF (ln_bdy) THEN 1234 ! Create a mask grid points in boundary rim 1235 IF (ld_bound_reject) THEN 1236 zbdymask(:,:) = 1.0_wp 1237 DO ji = 1, nb_bdy 1238 DO jj = 1, idx_bdy(ji)%nblen(1) 1239 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 1240 ENDDO 1241 ENDDO 1242 ENDIF 1243 1244 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 1245 ENDIF 1127 1246 1128 1247 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 1129 1248 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 1130 1249 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 1131 IF ( .NOT.( ln_zps .OR. ln_zco ) ) THEN 1132 ! Need to know the bathy depth for each observation for sco 1133 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 1250 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 1134 1251 & zgdepw ) 1135 ENDIF1136 1252 1137 1253 DO jobs = 1, kprofno 1138 1254 1139 1255 ! Skip bad profiles 1140 IF ( kpobsqc(jobs) >= 10) CYCLE1256 IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE 1141 1257 1142 1258 ! Check if this observation is on a grid point … … 1149 1265 IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 1150 1266 & .AND. & 1151 & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &1267 & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) & 1152 1268 & ) THEN 1153 1269 lgridobs = .TRUE. … … 1176 1292 & .OR. ( pobsdep(jobsp) < 0.0 ) & 1177 1293 & .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN 1178 kobsqc(jobsp) = kobsqc(jobsp) + 111294 kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 1179 1295 kosdobs = kosdobs + 1 1180 1296 CYCLE 1181 1297 ENDIF 1182 1298 1183 ! To check if an observations falls within land there are two cases: 1184 ! 1: z-coordibnates, where the check uses the mask 1185 ! 2: terrain following (eg s-coordinates), 1186 ! where we use the depth of the bottom cell to mask observations 1299 ! To check if an observations falls within land: 1187 1300 1188 IF( ln_zps .OR. ln_zco ) THEN !(CASE 1) 1301 ! Flag if the observation is deeper than the bathymetry 1302 ! Or if it is within the mask 1303 IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 1304 & .OR. & 1305 & ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1306 & == 0.0_wp) ) THEN 1307 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1308 klanobs = klanobs + 1 1309 CYCLE 1310 ENDIF 1189 1311 1190 ! Flag if the observation falls with a model land cell 1191 IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1192 & == 0.0_wp ) THEN 1193 kobsqc(jobsp) = kobsqc(jobsp) + 12 1194 klanobs = klanobs + 1 1195 CYCLE 1196 ENDIF 1197 1198 ! Flag if the observation is close to land 1199 IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 1200 & 0.0_wp) THEN 1201 knlaobs = knlaobs + 1 1202 IF (ld_nea) THEN 1203 kobsqc(jobsp) = kobsqc(jobsp) + 14 1204 ENDIF 1205 ENDIF 1206 1207 ELSE ! Case 2 1208 1209 ! Flag if the observation is deeper than the bathymetry 1210 ! Or if it is within the mask 1211 IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 1212 & .OR. & 1213 & ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1214 & == 0.0_wp) ) THEN 1215 kobsqc(jobsp) = kobsqc(jobsp) + 12 1216 klanobs = klanobs + 1 1217 CYCLE 1218 ENDIF 1219 1220 ! Flag if the observation is close to land 1221 IF ( ll_next_to_land ) THEN 1222 knlaobs = knlaobs + 1 1223 IF (ld_nea) THEN 1224 kobsqc(jobsp) = kobsqc(jobsp) + 14 1225 ENDIF 1226 ENDIF 1312 ! Flag if the observation is close to land 1313 IF ( ll_next_to_land ) THEN 1314 knlaobs = knlaobs + 1 1315 IF (ld_nea) THEN 1316 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1317 ENDIF 1227 1318 ENDIF 1228 1319 … … 1232 1323 IF (lgridobs) THEN 1233 1324 IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN 1234 kobsqc(jobsp) = kobsqc(jobsp) + 121325 kobsqc(jobsp) = IBSET(kobsqc(jobs),10) 1235 1326 klanobs = klanobs + 1 1236 1327 CYCLE … … 1250 1341 ENDIF 1251 1342 1343 IF (ln_bdy) THEN 1344 ! Flag if the observation falls close to the boundary rim 1345 IF (ld_bound_reject) THEN 1346 IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 1347 kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 1348 kbdyobs = kbdyobs + 1 1349 CYCLE 1350 ENDIF 1351 ! for observations on the grid... 1352 IF (lgridobs) THEN 1353 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 1354 kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 1355 kbdyobs = kbdyobs + 1 1356 CYCLE 1357 ENDIF 1358 ENDIF 1359 ENDIF 1360 ENDIF 1361 1252 1362 END DO 1253 1363 END DO … … 1255 1365 END SUBROUTINE obs_coo_spc_3d 1256 1366 1257 SUBROUTINE obs_pro_rej( profdata )1367 SUBROUTINE obs_pro_rej( profdata, kqc_cutoff ) 1258 1368 !!---------------------------------------------------------------------- 1259 1369 !! *** ROUTINE obs_pro_rej *** … … 1273 1383 !! * Arguments 1274 1384 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data 1385 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1386 1275 1387 !! * Local declarations 1276 1388 INTEGER :: jprof … … 1282 1394 DO jprof = 1, profdata%nprof 1283 1395 1284 IF ( profdata%nqc(jprof) > 10) THEN1396 IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN 1285 1397 1286 1398 DO jvar = 1, profdata%nvar … … 1290 1402 1291 1403 profdata%var(jvar)%nvqc(jobs) = & 1292 & profdata%var(jvar)%nvqc(jobs) + 261404 & IBSET(profdata%var(jvar)%nvqc(jobs),14) 1293 1405 1294 1406 END DO … … 1302 1414 END SUBROUTINE obs_pro_rej 1303 1415 1304 SUBROUTINE obs_uv_rej( profdata, knumu, knumv )1416 SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) 1305 1417 !!---------------------------------------------------------------------- 1306 1418 !! *** ROUTINE obs_uv_rej *** … … 1322 1434 INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected 1323 1435 INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected 1436 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1437 1324 1438 !! * Local declarations 1325 1439 INTEGER :: jprof … … 1341 1455 DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1) 1342 1456 1343 IF ( ( profdata%var(1)%nvqc(jobs) > 10) .AND. &1344 & ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN1345 profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 421457 IF ( ( profdata%var(1)%nvqc(jobs) > kqc_cutoff ) .AND. & 1458 & ( profdata%var(2)%nvqc(jobs) <= kqc_cutoff) ) THEN 1459 profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 1346 1460 knumv = knumv + 1 1347 1461 ENDIF 1348 IF ( ( profdata%var(2)%nvqc(jobs) > 10) .AND. &1349 & ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN1350 profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 421462 IF ( ( profdata%var(2)%nvqc(jobs) > kqc_cutoff ) .AND. & 1463 & ( profdata%var(1)%nvqc(jobs) <= kqc_cutoff) ) THEN 1464 profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 1351 1465 knumu = knumu + 1 1352 1466 ENDIF -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90
r6140 r8992 104 104 ! Bookkeeping arrays with sizes equal to number of variables 105 105 106 CHARACTER(len= 6), POINTER, DIMENSION(:) :: &106 CHARACTER(len=8), POINTER, DIMENSION(:) :: & 107 107 & cvars !: Variable names 108 108 -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
r6140 r8992 87 87 CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 88 88 CHARACTER(len=8) :: clrefdate 89 CHARACTER(len= 6), DIMENSION(:), ALLOCATABLE :: clvars89 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 90 90 INTEGER :: jvar 91 91 INTEGER :: ji … … 307 307 inowin = 0 308 308 DO ji = 1, inpfiles(jj)%nobs 309 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE310 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &311 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE309 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 310 IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 311 & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 312 312 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 313 313 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 325 325 inowin = 0 326 326 DO ji = 1, inpfiles(jj)%nobs 327 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE328 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &329 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE327 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 328 IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 329 & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 330 330 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 331 331 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 351 351 inowin = 0 352 352 DO ji = 1, inpfiles(jj)%nobs 353 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE354 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &355 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE353 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 354 IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 355 & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 356 356 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 357 357 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 373 373 374 374 DO ji = 1, inpfiles(jj)%nobs 375 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE376 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &377 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE375 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 376 IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 377 & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 378 378 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 379 379 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 388 388 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 389 389 & CYCLE 390 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &391 & ( inpfiles(jj)%idqc(ij,ji) <= 2) ) THEN390 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 391 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 392 392 ivar1t0 = ivar1t0 + 1 393 393 ENDIF … … 398 398 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 399 399 & CYCLE 400 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &401 & ( inpfiles(jj)%idqc(ij,ji) <= 2) ) THEN400 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 401 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 402 402 ivar2t0 = ivar2t0 + 1 403 403 ENDIF … … 407 407 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 408 408 & CYCLE 409 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &410 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) .AND. &411 & 412 & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &413 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) .AND. &409 IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 410 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 411 & ldvar1 ) .OR. & 412 & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 413 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 414 414 & ldvar2 ) ) THEN 415 415 ip3dt = ip3dt + 1 … … 437 437 DO jj = 1, inobf 438 438 DO ji = 1, inpfiles(jj)%nobs 439 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE440 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &441 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE439 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 440 IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 441 & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 442 442 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 443 443 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 452 452 DO jj = 1, inobf 453 453 DO ji = 1, inpfiles(jj)%nobs 454 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE455 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &456 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE454 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 455 IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 456 & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 457 457 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 458 458 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 501 501 ji = iprofidx(iindx(jk)) 502 502 503 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE504 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &505 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE503 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 504 IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 505 & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 506 506 507 507 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & … … 518 518 IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 519 519 520 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 521 & ( inpfiles(jj)%ivqc(ji,2) > 2 ) ) CYCLE 520 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 521 IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 522 & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 522 523 523 524 loop_prof : DO ij = 1, inpfiles(jj)%nlev … … 526 527 & CYCLE 527 528 528 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &529 & ( inpfiles(jj)%idqc(ij,ji) <= 2) ) THEN529 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 530 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 530 531 531 532 llvalprof = .TRUE. … … 534 535 ENDIF 535 536 536 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &537 & ( inpfiles(jj)%idqc(ij,ji) <= 2) ) THEN537 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 538 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 538 539 539 540 llvalprof = .TRUE. … … 615 616 IF (ldsatt) THEN 616 617 617 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &618 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) .AND. &619 & 620 & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &621 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) .AND. &622 & 618 IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 619 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 620 & ldvar1 ) .OR. & 621 & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 622 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 623 & ldvar2 ) ) THEN 623 624 ip3dt = ip3dt + 1 624 625 ELSE … … 628 629 ENDIF 629 630 630 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &631 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) .AND. &632 &ldvar1 ) .OR. ldsatt ) THEN631 IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 632 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 633 & ldvar1 ) .OR. ldsatt ) THEN 633 634 634 635 IF (ldsatt) THEN … … 661 662 662 663 ! Profile var1 value 663 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &664 & ( inpfiles(jj)%idqc(ij,ji) <= 2) ) THEN664 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 665 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 665 666 profdata%var(1)%vobs(ivar1t) = & 666 667 & inpfiles(jj)%pob(ij,ji,1) … … 692 693 ENDIF 693 694 694 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &695 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ).AND. &695 IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 696 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 696 697 & ldvar2 ) .OR. ldsatt ) THEN 697 698 … … 725 726 726 727 ! Profile var2 value 727 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &728 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) THEN728 IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) ) .AND. & 729 & ( .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) ) THEN 729 730 profdata%var(2)%vobs(ivar2t) = & 730 731 & inpfiles(jj)%pob(ij,ji,2) -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90
r6140 r8992 77 77 CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf' 78 78 CHARACTER(len=8) :: clrefdate 79 CHARACTER(len= 6), DIMENSION(:), ALLOCATABLE :: clvars79 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 80 80 INTEGER :: ji 81 81 INTEGER :: jj … … 172 172 173 173 !------------------------------------------------------------------ 174 ! Read the profile file into inpfiles174 ! Read the surface file into inpfiles 175 175 !------------------------------------------------------------------ 176 176 CALL init_obfbdata( inpfiles(jj) ) … … 296 296 ENDIF 297 297 llvalprof = .FALSE. 298 IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. & 299 & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN 298 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN 300 299 iobs = iobs + 1 301 300 ENDIF … … 370 369 ! Set observation information 371 370 372 IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. & 373 & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN 371 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN 374 372 375 373 iobs = iobs + 1 -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90
r6140 r8992 154 154 155 155 ! mark any masked data with a QC flag 156 IF( zobsmask(1) == 0 ) sladata%nqc(jobs) = 11156 IF( zobsmask(1) == 0 ) sladata%nqc(jobs) = IBSET(sladata%nqc(jobs),15) 157 157 158 158 END DO -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_sstbias.F90
r6140 r8992 1 1 MODULE obs_sstbias 2 2 !!====================================================================== 3 !! *** MODULE obs_ readsstbias ***4 !! Observation diagnostics: Read the bias for S LAdata3 !! *** MODULE obs_sstbias *** 4 !! Observation diagnostics: Read the bias for SST data 5 5 !!====================================================================== 6 6 !!---------------------------------------------------------------------- 7 !! obs_ rea_sstbias : Driver for reading altimeterbias7 !! obs_app_sstbias : Driver for reading and applying the SST bias 8 8 !!---------------------------------------------------------------------- 9 9 !! * Modules used … … 139 139 cl_bias_files(jtype) ) 140 140 ! Get the SST bias data 141 CALL iom_get( numsstbias, jpdom_data, ' sst', z_sstbias_2d(:,:), 1 )141 CALL iom_get( numsstbias, jpdom_data, 'tn', z_sstbias_2d(:,:), 1 ) 142 142 z_sstbias(:,:,jtype) = z_sstbias_2d(:,:) 143 143 ! Close the file … … 190 190 igrdj_tmp(:,:,jt) = igrdj(:,:,jobs) 191 191 zglam_tmp(:,:,jt) = zglam(:,:,jobs) 192 zgphi_tmp(:,:,jt) = zgphi(:,:,jobs)193 192 zgphi_tmp(:,:,jt) = zgphi(:,:,jobs) 194 193 zmask_tmp(:,:,jt) = zmask(:,:,jobs) -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90
r6140 r8992 50 50 INTEGER :: npj 51 51 INTEGER :: nsurfup !: Observation counter used in obs_oper 52 INTEGER :: nrec !: Number of surface observation records in window 52 53 53 54 ! Arrays with size equal to the number of surface observations … … 56 57 & mi, & !: i-th grid coord. for interpolating to surface observation 57 58 & mj, & !: j-th grid coord. for interpolating to surface observation 59 & mt, & !: time record number for gridded data 58 60 & nsidx,& !: Surface observation number 59 61 & nsfil,& !: Surface observation number in file … … 67 69 & ntyp !: Type of surface observation product 68 70 69 CHARACTER(len= 6), POINTER, DIMENSION(:) :: &71 CHARACTER(len=8), POINTER, DIMENSION(:) :: & 70 72 & cvars !: Variable names 71 73 … … 93 95 & nsstpmpp !: Global number of surface observations per time step 94 96 97 ! Arrays with size equal to the number of observation records in the window 98 INTEGER, POINTER, DIMENSION(:) :: & 99 & mrecstp ! Time step of the records 100 95 101 ! Arrays used to store source indices when 96 102 ! compressing obs_surf derived types … … 100 106 INTEGER, POINTER, DIMENSION(:) :: & 101 107 & nsind !: Source indices of surface data in compressed data 108 109 ! Is this a gridded product? 110 111 LOGICAL :: lgrid 102 112 103 113 END TYPE obs_surf … … 160 170 & surf%mi(ksurf), & 161 171 & surf%mj(ksurf), & 172 & surf%mt(ksurf), & 162 173 & surf%nsidx(ksurf), & 163 174 & surf%nsfil(ksurf), & … … 176 187 & ) 177 188 189 surf%mt(:) = -1 190 178 191 179 192 ! Allocate arrays of number of surface data size * number of variables … … 190 203 & ) 191 204 205 surf%rext(:,:) = 0.0_wp 206 192 207 ! Allocate arrays of number of time step size 193 208 … … 217 232 218 233 surf%nsurfup = 0 234 235 ! Not gridded by default 236 237 surf%lgrid = .FALSE. 219 238 220 239 END SUBROUTINE obs_surf_alloc … … 242 261 & surf%mi, & 243 262 & surf%mj, & 263 & surf%mt, & 244 264 & surf%nsidx, & 245 265 & surf%nsfil, & … … 370 390 newsurf%mi(insurf) = surf%mi(ji) 371 391 newsurf%mj(insurf) = surf%mj(ji) 392 newsurf%mt(insurf) = surf%mt(ji) 372 393 newsurf%nsidx(insurf) = surf%nsidx(ji) 373 394 newsurf%nsfil(insurf) = surf%nsfil(ji) … … 414 435 newsurf%nstp = surf%nstp 415 436 newsurf%cvars(:) = surf%cvars(:) 437 438 ! Set gridded stuff 439 440 newsurf%mt(insurf) = surf%mt(ji) 416 441 417 442 ! Deallocate temporary data … … 454 479 oldsurf%mi(jj) = surf%mi(ji) 455 480 oldsurf%mj(jj) = surf%mj(ji) 481 oldsurf%mt(jj) = surf%mt(ji) 456 482 oldsurf%nsidx(jj) = surf%nsidx(ji) 457 483 oldsurf%nsfil(jj) = surf%nsfil(ji) -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
r6140 r8992 8 8 !! obs_wri_prof : Write profile observations in feedback format 9 9 !! obs_wri_surf : Write surface observations in feedback format 10 !! obs_wri_stats : Print basic statistics on the data being written out10 !! obs_wri_stats : Print basic statistics on the data being written out 11 11 !!---------------------------------------------------------------------- 12 12 … … 83 83 TYPE(obfbdata) :: fbdata 84 84 CHARACTER(LEN=40) :: clfname 85 CHARACTER(LEN= 6) :: clfiletype85 CHARACTER(LEN=10) :: clfiletype 86 86 INTEGER :: ilevel 87 87 INTEGER :: jvar … … 196 196 fbdata%ivqc(jo,:) = profdata%ivqc(jo,:) 197 197 fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:) 198 IF ( profdata%nqc(jo) > 10) THEN199 fbdata%ioqc(jo) = 4198 IF ( profdata%nqc(jo) > 255 ) THEN 199 fbdata%ioqc(jo) = IBSET(profdata%nqc(jo),2) 200 200 fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo) 201 fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10201 fbdata%ioqcf(2,jo) = profdata%nqc(jo) 202 202 ELSE 203 203 fbdata%ioqc(jo) = profdata%nqc(jo) … … 236 236 fbdata%idqc(ik,jo) = profdata%var(jvar)%idqc(jk) 237 237 fbdata%idqcf(:,ik,jo) = profdata%var(jvar)%idqcf(:,jk) 238 IF ( profdata%var(jvar)%nvqc(jk) > 10) THEN239 fbdata%ivlqc(ik,jo,jvar) = 4238 IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN 239 fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2) 240 240 fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk) 241 fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10 241 !$AGRIF_DO_NOT_TREAT 242 fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000 0000 1111 1111') 243 !$AGRIF_END_DO_NOT_TREAT 242 244 ELSE 243 245 fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) … … 320 322 TYPE(obfbdata) :: fbdata 321 323 CHARACTER(LEN=40) :: clfname ! netCDF filename 322 CHARACTER(LEN= 6):: clfiletype324 CHARACTER(LEN=10) :: clfiletype 323 325 CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf' 324 326 INTEGER :: jo … … 395 397 END DO 396 398 397 CASE('ICECON ')399 CASE('ICECONC') 398 400 399 401 CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & … … 418 420 END DO 419 421 422 CASE('SSS') 423 424 CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 425 & 1 + iadd, iext, .TRUE. ) 426 427 clfiletype = 'sssfb' 428 fbdata%cname(1) = surfdata%cvars(1) 429 fbdata%coblong(1) = 'Sea surface salinity' 430 fbdata%cobunit(1) = 'psu' 431 DO je = 1, iext 432 fbdata%cextname(je) = pext%cdname(je) 433 fbdata%cextlong(je) = pext%cdlong(je,1) 434 fbdata%cextunit(je) = pext%cdunit(je,1) 435 END DO 436 fbdata%caddlong(1,1) = 'Model interpolated SSS' 437 fbdata%caddunit(1,1) = 'psu' 438 fbdata%cgrid(1) = 'T' 439 DO ja = 1, iadd 440 fbdata%caddname(1+ja) = padd%cdname(ja) 441 fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 442 fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 443 END DO 444 445 CASE DEFAULT 446 447 CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' ) 448 420 449 END SELECT 421 450 … … 439 468 fbdata%ivqc(jo,:) = 0 440 469 fbdata%ivqcf(:,jo,:) = 0 441 IF ( surfdata%nqc(jo) > 10) THEN470 IF ( surfdata%nqc(jo) > 255 ) THEN 442 471 fbdata%ioqc(jo) = 4 443 472 fbdata%ioqcf(1,jo) = 0 444 fbdata%ioqcf(2,jo) = surfdata%nqc(jo) - 10 473 !$AGRIF_DO_NOT_TREAT 474 fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111') 475 !$AGRIF_END_DO_NOT_TREAT 445 476 ELSE 446 477 fbdata%ioqc(jo) = surfdata%nqc(jo) … … 474 505 fbdata%idqc(1,jo) = 0 475 506 fbdata%idqcf(:,1,jo) = 0 476 IF ( surfdata%nqc(jo) > 10) THEN507 IF ( surfdata%nqc(jo) > 255 ) THEN 477 508 fbdata%ivqc(jo,1) = 4 478 509 fbdata%ivlqc(1,jo,1) = 4 479 510 fbdata%ivlqcf(1,1,jo,1) = 0 480 fbdata%ivlqcf(2,1,jo,1) = surfdata%nqc(jo) - 10 511 !$AGRIF_DO_NOT_TREAT 512 fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111') 513 !$AGRIF_END_DO_NOT_TREAT 481 514 ELSE 482 515 fbdata%ivqc(jo,1) = surfdata%nqc(jo) -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_h2d.h90
r2474 r8992 1240 1240 & zdum, & 1241 1241 & zaamax 1242 1242 1243 imax = -1 1243 1244 ! Main computation 1244 1245 pflt = 1.0_wp -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90
r7851 r8992 66 66 INTEGER :: nsnd ! total number of fields sent 67 67 INTEGER :: ncplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 68 INTEGER, PUBLIC, PARAMETER :: nmaxfld= 55! Maximum number of coupling fields68 INTEGER, PUBLIC, PARAMETER :: nmaxfld=60 ! Maximum number of coupling fields 69 69 INTEGER, PUBLIC, PARAMETER :: nmaxcat=5 ! Maximum number of coupling fields 70 70 INTEGER, PUBLIC, PARAMETER :: nmaxcpl=5 ! Maximum number of coupling fields -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r7788 r8992 65 65 LOGICAL , PUBLIC :: ln_sdw !: true if 3d stokes drift from wave model 66 66 LOGICAL , PUBLIC :: ln_tauoc !: true if normalized stress from wave is used 67 LOGICAL , PUBLIC :: ln_tauw !: true if ocean stress components from wave is used 67 68 LOGICAL , PUBLIC :: ln_stcor !: true if Stokes-Coriolis term is used 69 ! 70 INTEGER , PUBLIC :: nn_sdrift ! type of parameterization to calculate vertical Stokes drift 68 71 ! 69 72 LOGICAL , PUBLIC :: ln_icebergs !: Icebergs … … 79 82 INTEGER , PUBLIC, PARAMETER :: jp_none = 5 !: for OPA when doing coupling via SAS module 80 83 84 !!---------------------------------------------------------------------- 85 !! Stokes drift parametrization definition 86 !!---------------------------------------------------------------------- 87 INTEGER , PUBLIC, PARAMETER :: jp_breivik = 0 ! Breivik 2015: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 88 INTEGER , PUBLIC, PARAMETER :: jp_phillips = 1 ! Phillips: v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))] 89 INTEGER , PUBLIC, PARAMETER :: jp_peakfr = 2 ! Phillips using the peak wave number read from wave model instead of the inverse depth scale 90 81 91 !!---------------------------------------------------------------------- 82 92 !! component definition -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
r8329 r8992 113 113 INTEGER, PARAMETER :: jpr_wper = 48 ! Mean wave period 114 114 INTEGER, PARAMETER :: jpr_wnum = 49 ! Mean wavenumber 115 INTEGER, PARAMETER :: jpr_ wstrf= 50 ! Stress fraction adsorbed by waves115 INTEGER, PARAMETER :: jpr_tauoc = 50 ! Stress fraction adsorbed by waves 116 116 INTEGER, PARAMETER :: jpr_wdrag = 51 ! Neutral surface drag coefficient 117 117 INTEGER, PARAMETER :: jpr_isf = 52 118 118 INTEGER, PARAMETER :: jpr_icb = 53 119 120 INTEGER, PARAMETER :: jprcv = 53 ! total number of fields received 119 INTEGER, PARAMETER :: jpr_wfreq = 54 ! Wave peak frequency 120 INTEGER, PARAMETER :: jpr_tauwx = 55 ! x component of the ocean stress from waves 121 INTEGER, PARAMETER :: jpr_tauwy = 56 ! y component of the ocean stress from waves 122 123 INTEGER, PARAMETER :: jprcv = 56 ! total number of fields received 121 124 122 125 INTEGER, PARAMETER :: jps_fice = 1 ! ice fraction sent to the atmosphere … … 165 168 TYPE(FLD_C) :: sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2 166 169 ! ! Received from the atmosphere 167 TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf 170 TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_tauw, sn_rcv_dqnsdt, sn_rcv_qsr, & 171 sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf 168 172 TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf 169 173 ! Send to waves 170 174 TYPE(FLD_C) :: sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 171 175 ! Received from waves 172 TYPE(FLD_C) :: sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_wstrf,sn_rcv_wdrag 176 TYPE(FLD_C) :: sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_tauoc,sn_rcv_wdrag, & 177 sn_rcv_wfreq 173 178 ! ! Other namelist parameters 174 179 INTEGER :: nn_cplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data … … 242 247 & sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau , sn_rcv_dqnsdt, sn_rcv_qsr, & 243 248 & sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev , sn_rcv_hsig , sn_rcv_phioc , & 244 & sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper , sn_rcv_wnum , sn_rcv_ wstrf, &249 & sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper , sn_rcv_wnum , sn_rcv_tauoc , & 245 250 & sn_rcv_wdrag, sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf , sn_rcv_cal , & 246 251 & sn_rcv_iceflx,sn_rcv_co2 , nn_cplmodel , ln_usecplmask, sn_rcv_mslp , & 247 & sn_rcv_icb , sn_rcv_isf 252 & sn_rcv_icb , sn_rcv_isf , sn_rcv_wfreq , sn_rcv_tauw 248 253 249 254 !!--------------------------------------------------------------------- … … 295 300 WRITE(numout,*)' Mean wave period = ', TRIM(sn_rcv_wper%cldes ), ' (', TRIM(sn_rcv_wper%clcat ), ')' 296 301 WRITE(numout,*)' Mean wave number = ', TRIM(sn_rcv_wnum%cldes ), ' (', TRIM(sn_rcv_wnum%clcat ), ')' 297 WRITE(numout,*)' Stress frac adsorbed by waves = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')' 302 WRITE(numout,*)' Wave peak frequency = ', TRIM(sn_rcv_wfreq%cldes ), ' (', TRIM(sn_rcv_wfreq%clcat ), ')' 303 WRITE(numout,*)' Stress frac adsorbed by waves = ', TRIM(sn_rcv_tauoc%cldes ), ' (', TRIM(sn_rcv_tauoc%clcat ), ')' 304 WRITE(numout,*)' Stress components by waves = ', TRIM(sn_rcv_tauw%cldes ), ' (', TRIM(sn_rcv_tauw%clcat ), ')' 298 305 WRITE(numout,*)' Neutral surf drag coefficient = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 299 306 WRITE(numout,*)' sent fields (multiple ice categories)' … … 578 585 cpl_wper = .TRUE. 579 586 ENDIF 587 srcv(jpr_wfreq)%clname = 'O_WFreq' ! wave peak frequency 588 IF( TRIM(sn_rcv_wfreq%cldes ) == 'coupled' ) THEN 589 srcv(jpr_wfreq)%laction = .TRUE. 590 cpl_wfreq = .TRUE. 591 ENDIF 580 592 srcv(jpr_wnum)%clname = 'O_WNum' ! mean wave number 581 593 IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' ) THEN … … 583 595 cpl_wnum = .TRUE. 584 596 ENDIF 585 srcv(jpr_wstrf)%clname = 'O_WStrf' ! stress fraction adsorbed by the wave 586 IF( TRIM(sn_rcv_wstrf%cldes ) == 'coupled' ) THEN 587 srcv(jpr_wstrf)%laction = .TRUE. 588 cpl_wstrf = .TRUE. 597 srcv(jpr_tauoc)%clname = 'O_TauOce' ! stress fraction adsorbed by the wave 598 IF( TRIM(sn_rcv_tauoc%cldes ) == 'coupled' ) THEN 599 srcv(jpr_tauoc)%laction = .TRUE. 600 cpl_tauoc = .TRUE. 601 ENDIF 602 srcv(jpr_tauwx)%clname = 'O_Tauwx' ! ocean stress from wave in the x direction 603 srcv(jpr_tauwy)%clname = 'O_Tauwy' ! ocean stress from wave in the y direction 604 IF( TRIM(sn_rcv_tauw%cldes ) == 'coupled' ) THEN 605 srcv(jpr_tauwx)%laction = .TRUE. 606 srcv(jpr_tauwy)%laction = .TRUE. 607 cpl_tauw = .TRUE. 589 608 ENDIF 590 609 srcv(jpr_wdrag)%clname = 'O_WDrag' ! neutral surface drag coefficient … … 594 613 ENDIF 595 614 ! 615 IF( srcv(jpr_tauoc)%laction .AND. srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction ) & 616 CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', & 617 '(sn_rcv_tauoc=coupled and sn_rcv_tauw=coupled)' ) 618 ! 596 619 ! ! ------------------------------- ! 597 620 ! ! OPA-SAS coupling - rcv by opa ! … … 1165 1188 ! ! ========================= ! 1166 1189 IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1) 1190 ! 1191 ! ! ========================= ! 1192 ! ! Wave peak frequency ! 1193 ! ! ========================= ! 1194 IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1) 1167 1195 ! 1168 1196 ! ! ========================= ! … … 1173 1201 ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode 1174 1202 IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction & 1175 .OR. srcv(jpr_hsig)%laction ) THEN1203 .OR. srcv(jpr_hsig)%laction .OR. srcv(jpr_wfreq)%laction ) THEN 1176 1204 CALL sbc_stokes() 1177 1205 ENDIF … … 1180 1208 ! ! Stress adsorbed by waves ! 1181 1209 ! ! ========================= ! 1182 IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1) 1210 IF( srcv(jpr_tauoc)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_tauoc)%z3(:,:,1) 1211 1212 ! ! ========================= ! 1213 ! ! Stress component by waves ! 1214 ! ! ========================= ! 1215 IF( srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction .AND. ln_tauw ) THEN 1216 tauw_x(:,:) = frcv(jpr_tauwx)%z3(:,:,1) 1217 tauw_y(:,:) = frcv(jpr_tauwy)%z3(:,:,1) 1218 ENDIF 1183 1219 1184 1220 ! ! ========================= ! -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r8524 r8992 96 96 & ln_rnf , nn_fwb , ln_ssr , ln_isf , ln_apr_dyn , & 97 97 & ln_wave , ln_cdgw , ln_sdw , ln_tauoc , ln_stcor , & 98 & nn_lsm98 & ln_tauw , nn_lsm, nn_sdrift 99 99 !!---------------------------------------------------------------------- 100 100 ! … … 157 157 WRITE(numout,*) ' surface wave ln_wave = ', ln_wave 158 158 WRITE(numout,*) ' Stokes drift corr. to vert. velocity ln_sdw = ', ln_sdw 159 WRITE(numout,*) ' vertical parametrization nn_sdrift = ', nn_sdrift 159 160 WRITE(numout,*) ' wave modified ocean stress ln_tauoc = ', ln_tauoc 161 WRITE(numout,*) ' wave modified ocean stress component ln_tauw = ', ln_tauw 160 162 WRITE(numout,*) ' Stokes coriolis term ln_stcor = ', ln_stcor 161 163 WRITE(numout,*) ' neutral drag coefficient (CORE, MFS) ln_cdgw = ', ln_cdgw 162 164 ENDIF 165 ! 166 IF( ln_sdw ) THEN 167 IF( .NOT.(nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakfr) ) & 168 CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' ) 169 ENDIF 170 IF( ln_tauoc .AND. ln_tauw ) & 171 CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', & 172 '(ln_tauoc=.true. and ln_tauw=.true.)' ) 173 IF( ln_tauoc ) & 174 CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauoc=.true.)' ) 175 IF( ln_tauw ) & 176 CALL ctl_warn( 'The wave modified ocean stress components are used (ln_tauw=.true.) ', & 177 'This will override any other specification of the ocean stress' ) 163 178 ! 164 179 IF( .NOT.ln_usr ) THEN ! the model calendar needs some specificities (except in user defined case) … … 410 425 IF( ll_opa ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice ) ! OPA-SAS coupling: OPA receiving fields from SAS 411 426 END SELECT 412 IF ( ln_wave .AND. ln_tauoc) THEN ! Wave stress subctracted 413 utau(:,:) = utau(:,:)*tauoc_wave(:,:) 414 vtau(:,:) = vtau(:,:)*tauoc_wave(:,:) 415 taum(:,:) = taum(:,:)*tauoc_wave(:,:) 416 ! 417 SELECT CASE( nsbc ) 418 CASE( 0,1,2,3,5,-1 ) ; 419 IF(lwp .AND. kt == nit000 ) WRITE(numout,*) 'WARNING: You are subtracting the wave stress to the ocean. & 420 & If not requested select ln_tauoc=.false' 421 END SELECT 422 ! 423 END IF 427 ! 424 428 IF( ln_mixcpl ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice ) ! forced-coupled mixed formulation after forcing 425 429 ! 430 IF ( ln_wave .AND. (ln_tauoc .OR. ln_tauw) ) CALL sbc_wstress( ) ! Wind stress provided by waves 426 431 ! 427 432 ! !== Misc. Options ==! -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
r7968 r8992 39 39 ! !!* namsbc_rnf namelist * 40 40 CHARACTER(len=100) :: cn_dir !: Root directory for location of rnf files 41 LOGICAL 41 LOGICAL , PUBLIC :: ln_rnf_depth !: depth river runoffs attribute specified in a file 42 42 LOGICAL :: ln_rnf_depth_ini !: depth river runoffs computed at the initialisation 43 43 REAL(wp) :: rn_rnf_max !: maximum value of the runoff climatologie (ln_rnf_depth_ini =T) … … 122 122 IF( .NOT. l_rnfcpl ) rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) ! updated runoff value at time step kt 123 123 ! 124 ! ! set temperature & salinity content of runoffs124 ! ! set temperature & salinity content of runoffs 125 125 IF( ln_rnf_tem ) THEN ! use runoffs temperature data 126 126 rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 … … 133 133 END WHERE 134 134 ELSE ! use SST as runoffs temperature 135 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 135 !CEOD River is fresh water so must at least be 0 unless we consider ice 136 rnf_tsc(:,:,jp_tem) = MAX(sst_m(:,:),0.0_wp) * rnf(:,:) * r1_rau0 136 137 ENDIF 137 138 ! ! use runoffs salinity data -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r7864 r8992 33 33 34 34 PUBLIC sbc_stokes ! routine called in sbccpl 35 PUBLIC sbc_wstress ! routine called in sbcmod 35 36 PUBLIC sbc_wave ! routine called in sbcmod 36 37 PUBLIC sbc_wave_init ! routine called in sbcmod … … 42 43 LOGICAL, PUBLIC :: cpl_sdrfty = .FALSE. 43 44 LOGICAL, PUBLIC :: cpl_wper = .FALSE. 45 LOGICAL, PUBLIC :: cpl_wfreq = .FALSE. 44 46 LOGICAL, PUBLIC :: cpl_wnum = .FALSE. 45 LOGICAL, PUBLIC :: cpl_wstrf = .FALSE. 47 LOGICAL, PUBLIC :: cpl_tauoc = .FALSE. 48 LOGICAL, PUBLIC :: cpl_tauw = .FALSE. 46 49 LOGICAL, PUBLIC :: cpl_wdrag = .FALSE. 47 50 … … 51 54 INTEGER :: jp_hsw ! index of significant wave hight (m) at T-point 52 55 INTEGER :: jp_wmp ! index of mean wave period (s) at T-point 56 INTEGER :: jp_wfr ! index of wave peak frequency (1/s) at T-point 53 57 54 58 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_cd ! structure of input fields (file informations, fields read) Drag Coefficient … … 56 60 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_wn ! structure of input fields (file informations, fields read) wave number for Qiao 57 61 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 62 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tauw ! structure of input fields (file informations, fields read) ocean stress components from wave model 63 58 64 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: cdn_wave !: 59 65 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hsw, wmp, wnum !: 66 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wfreq !: 60 67 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tauoc_wave !: 68 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tauw_x, tauw_y !: 61 69 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tsd2d !: 62 70 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: div_sd !: barotropic stokes drift divergence … … 96 104 CALL wrk_alloc( jpi,jpj, zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 97 105 ! 98 ! 99 zfac = 2.0_wp * rpi / 16.0_wp 100 DO jj = 1, jpj ! exp. wave number at t-point (Eq. (19) in Breivick et al. (2014) ) 101 DO ji = 1, jpi 106 ! select parameterization for the calculation of vertical Stokes drift 107 ! exp. wave number at t-point 108 IF( nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips ) THEN ! (Eq. (19) in Breivick et al. (2014) ) 109 zfac = 2.0_wp * rpi / 16.0_wp 110 DO jj = 1, jpj 111 DO ji = 1, jpi 102 112 ! Stokes drift velocity estimated from Hs and Tmean 103 ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj) 113 ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp ) 104 114 ! Stokes surface speed 105 zsp0 = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj) ) 106 tsd2d(ji,jj) = zsp0 115 tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj)) 107 116 ! Wavenumber scale 108 zk_t(ji,jj) = ABS( zsp0 ) / MAX( ABS( 5.97_wp*ztransp ) , 0.0000001_wp ) 109 END DO 110 END DO 111 DO jj = 1, jpjm1 ! exp. wave number & Stokes drift velocity at u- & v-points 112 DO ji = 1, jpim1 113 zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 114 zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 115 ! 116 zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 117 zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 118 END DO 119 END DO 117 zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp ) 118 END DO 119 END DO 120 DO jj = 1, jpjm1 ! exp. wave number & Stokes drift velocity at u- & v-points 121 DO ji = 1, jpim1 122 zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 123 zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 124 ! 125 zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 126 zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 127 END DO 128 END DO 129 ELSE IF( nn_sdrift==jp_peakfr ) THEN ! peak wave number calculated from the peak frequency received by the wave model 130 DO jj = 1, jpjm1 131 DO ji = 1, jpim1 132 zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav 133 zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav 134 ! 135 zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 136 zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 137 END DO 138 END DO 139 ENDIF 120 140 ! 121 141 ! !== horizontal Stokes Drift 3D velocity ==! 122 DO jk = 1, jpkm1 123 DO jj = 2, jpjm1 124 DO ji = 2, jpim1 125 zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 126 zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 127 ! 128 zkh_u = zk_u(ji,jj) * zdep_u ! k * depth 129 zkh_v = zk_v(ji,jj) * zdep_v 130 ! ! Depth attenuation 131 zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u ) 132 zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v ) 133 ! 134 usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 135 vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 136 END DO 137 END DO 138 END DO 142 IF( nn_sdrift==jp_breivik ) THEN 143 DO jk = 1, jpkm1 144 DO jj = 2, jpjm1 145 DO ji = 2, jpim1 146 zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 147 zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 148 ! 149 zkh_u = zk_u(ji,jj) * zdep_u ! k * depth 150 zkh_v = zk_v(ji,jj) * zdep_v 151 ! ! Depth attenuation 152 zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u ) 153 zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v ) 154 ! 155 usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 156 vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 157 END DO 158 END DO 159 END DO 160 ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakfr ) THEN 161 DO jk = 1, jpkm1 162 DO jj = 2, jpjm1 163 DO ji = 2, jpim1 164 zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 165 zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 166 ! 167 zkh_u = zk_u(ji,jj) * zdep_u ! k * depth 168 zkh_v = zk_v(ji,jj) * zdep_v 169 ! ! Depth attenuation 170 zda_u = EXP( -2.0_wp*zkh_u ) - SQRT(2.0_wp*rpi*zkh_u) * ERFC(SQRT(2.0_wp*zkh_u)) 171 zda_v = EXP( -2.0_wp*zkh_v ) - SQRT(2.0_wp*rpi*zkh_v) * ERFC(SQRT(2.0_wp*zkh_v)) 172 ! 173 usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 174 vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 175 END DO 176 END DO 177 END DO 178 ENDIF 179 139 180 CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. ) 140 181 ! … … 189 230 190 231 232 SUBROUTINE sbc_wstress( ) 233 !!--------------------------------------------------------------------- 234 !! *** ROUTINE sbc_wstress *** 235 !! 236 !! ** Purpose : Updates the ocean momentum modified by waves 237 !! 238 !! ** Method : - Calculate u,v components of stress depending on stress 239 !! model 240 !! - Calculate the stress module 241 !! - The wind module is not modified by waves 242 !! ** action 243 !!--------------------------------------------------------------------- 244 INTEGER :: jj, ji ! dummy loop argument 245 ! 246 IF( ln_tauoc ) THEN 247 utau(:,:) = utau(:,:)*tauoc_wave(:,:) 248 vtau(:,:) = vtau(:,:)*tauoc_wave(:,:) 249 taum(:,:) = taum(:,:)*tauoc_wave(:,:) 250 ENDIF 251 ! 252 IF( ln_tauw ) THEN 253 DO jj = 1, jpjm1 254 DO ji = 1, jpim1 255 ! Stress components at u- & v-points 256 utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) ) 257 vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) ) 258 ! 259 ! Stress module at t points 260 taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) ) 261 END DO 262 END DO 263 CALL lbc_lnk_multi( utau(:,:), 'U', -1. , vtau(:,:), 'V', -1. , taum(:,:) , 'T', -1. ) 264 ENDIF 265 ! 266 END SUBROUTINE sbc_wstress 267 268 191 269 SUBROUTINE sbc_wave( kt ) 192 270 !!--------------------------------------------------------------------- … … 211 289 ENDIF 212 290 213 IF( ln_tauoc .AND. .NOT. cpl_ wstrf) THEN !== Wave induced stress ==!291 IF( ln_tauoc .AND. .NOT. cpl_tauoc ) THEN !== Wave induced stress ==! 214 292 CALL fld_read( kt, nn_fsbc, sf_tauoc ) ! read wave norm stress from external forcing 215 293 tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 294 ENDIF 295 296 IF( ln_tauw .AND. .NOT. cpl_tauw ) THEN !== Wave induced stress ==! 297 CALL fld_read( kt, nn_fsbc, sf_tauw ) ! read ocean stress components from external forcing (T grid) 298 tauw_x(:,:) = sf_tauw(1)%fnow(:,:,1) 299 tauw_y(:,:) = sf_tauw(2)%fnow(:,:,1) 216 300 ENDIF 217 301 … … 222 306 IF( jp_hsw > 0 ) hsw (:,:) = sf_sd(jp_hsw)%fnow(:,:,1) ! significant wave height 223 307 IF( jp_wmp > 0 ) wmp (:,:) = sf_sd(jp_wmp)%fnow(:,:,1) ! wave mean period 308 IF( jp_wfr > 0 ) wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1) ! Peak wave frequency 224 309 IF( jp_usd > 0 ) ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1) ! 2D zonal Stokes Drift at T point 225 310 IF( jp_vsd > 0 ) vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) ! 2D meridional Stokes Drift at T point … … 234 319 ! !== Computation of the 3d Stokes Drift ==! 235 320 ! 236 IF( jpfld == 4 ) CALL sbc_stokes() ! Calculate only if required fields are read 237 ! ! In coupled wave model-NEMO case the call is done after coupling 321 IF( ((nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) .AND. & 322 jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. & 323 (nn_sdrift==jp_peakfr .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) & 324 CALL sbc_stokes() ! Calculate only if required fields are read 325 ! ! In coupled wave model-NEMO case the call is done after coupling 238 326 ! 239 327 ENDIF … … 260 348 !! 261 349 CHARACTER(len=100) :: cn_dir ! Root directory for location of drag coefficient files 262 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read350 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i, slf_j ! array of namelist informations on the fields to read 263 351 TYPE(FLD_N) :: sn_cdg, sn_usd, sn_vsd, & 264 & sn_hsw, sn_wmp, sn_wnum, sn_tauoc ! informations about the fields to be read 265 ! 266 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc 352 & sn_hsw, sn_wmp, sn_wfr, sn_wnum, & 353 & sn_tauoc, sn_tauwx, sn_tauwy ! informations about the fields to be read 354 ! 355 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, & 356 sn_wnum, sn_tauoc, sn_tauwx, sn_tauwy 267 357 !!--------------------------------------------------------------------- 268 358 ! … … 289 379 290 380 IF( ln_tauoc ) THEN 291 IF( .NOT. cpl_ wstrf) THEN381 IF( .NOT. cpl_tauoc ) THEN 292 382 ALLOCATE( sf_tauoc(1), STAT=ierror ) !* allocate and fill sf_wave with sn_tauoc 293 383 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) … … 300 390 ENDIF 301 391 392 IF( ln_tauw ) THEN 393 IF( .NOT. cpl_tauw ) THEN 394 ALLOCATE( sf_tauw(2), STAT=ierror ) !* allocate and fill sf_wave with sn_tauwx/y 395 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauw structure' ) 396 ! 397 ALLOCATE( slf_j(2) ) 398 slf_j(1) = sn_tauwx 399 slf_j(2) = sn_tauwy 400 ALLOCATE( sf_tauw(1)%fnow(jpi,jpj,1) ) 401 ALLOCATE( sf_tauw(2)%fnow(jpi,jpj,1) ) 402 IF( slf_j(1)%ln_tint ) ALLOCATE( sf_tauw(1)%fdta(jpi,jpj,1,2) ) 403 IF( slf_j(2)%ln_tint ) ALLOCATE( sf_tauw(2)%fdta(jpi,jpj,1,2) ) 404 CALL fld_fill( sf_tauw, (/ slf_j /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 405 ENDIF 406 ALLOCATE( tauw_x(jpi,jpj) ) 407 ALLOCATE( tauw_y(jpi,jpj) ) 408 ENDIF 409 302 410 IF( ln_sdw ) THEN ! Find out how many fields have to be read from file if not coupled 303 411 jpfld=0 304 jp_usd=0 ; jp_vsd=0 ; jp_hsw=0 ; jp_wmp=0 412 jp_usd=0 ; jp_vsd=0 ; jp_hsw=0 ; jp_wmp=0 ; jp_wfr=0 305 413 IF( .NOT. cpl_sdrftx ) THEN 306 414 jpfld = jpfld + 1 … … 311 419 jp_vsd = jpfld 312 420 ENDIF 313 IF( .NOT. cpl_hsig ) THEN421 IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 314 422 jpfld = jpfld + 1 315 423 jp_hsw = jpfld 316 424 ENDIF 317 IF( .NOT. cpl_wper ) THEN425 IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 318 426 jpfld = jpfld + 1 319 427 jp_wmp = jpfld 428 ENDIF 429 IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakfr ) THEN 430 jpfld = jpfld + 1 431 jp_wfr = jpfld 320 432 ENDIF 321 433 … … 327 439 IF( jp_hsw > 0 ) slf_i(jp_hsw) = sn_hsw 328 440 IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 441 IF( jp_wfr > 0 ) slf_i(jp_wfr) = sn_wfr 442 329 443 ALLOCATE( sf_sd(jpfld), STAT=ierror ) !* allocate and fill sf_sd with stokes drift 330 444 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) … … 339 453 ALLOCATE( usd (jpi,jpj,jpk), vsd (jpi,jpj,jpk), wsd(jpi,jpj,jpk) ) 340 454 ALLOCATE( hsw (jpi,jpj) , wmp (jpi,jpj) ) 455 ALLOCATE( wfreq(jpi,jpj) ) 341 456 ALLOCATE( ut0sd(jpi,jpj) , vt0sd(jpi,jpj) ) 342 457 ALLOCATE( div_sd(jpi,jpj) ) -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
r8698 r8992 317 317 IF( jk == mikt(ji,jj) ) THEN ! first level 318 318 ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj) - emp(ji,jj) ) & 319 & - (rnf_b(ji,jj) - rnf(ji,jj) ) &320 319 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) 321 320 ztc_f = ztc_f - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 322 321 ENDIF 322 IF( ln_rnf_depth ) THEN 323 ! Rivers are not just at the surface must go down to nk_rnf(ji,jj) 324 IF( mikt(ji,jj) <=jk .and. jk <= nk_rnf(ji,jj) ) THEN 325 ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj) - rnf(ji,jj) ) ) & 326 & * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) 327 ENDIF 328 ELSE 329 IF( jk == mikt(ji,jj) ) THEN ! first level 330 ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj) - rnf(ji,jj) ) ) 331 ENDIF 332 ENDIF 333 323 334 ! 324 335 ! solar penetration (temperature only) -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
r7788 r8992 27 27 USE trd_oce ! trends: ocean variables 28 28 USE trdtra ! trends manager: tracers 29 #if defined key_asminc 30 USE asminc ! Assimilation increment 31 #endif 29 32 ! 30 33 USE in_out_manager ! I/O manager … … 34 37 USE wrk_nemo ! Memory Allocation 35 38 USE timing ! Timing 39 USE wet_dry, ONLY : ll_wd, rn_wdmin1, r_rn_wdmin1 ! Wetting and drying 36 40 37 41 IMPLICIT NONE … … 72 76 INTEGER, INTENT(in) :: kt ! ocean time-step index 73 77 ! 74 INTEGER :: ji, jj, jk, jn ! dummy loop indices75 INTEGER :: ikt, ikb ! local integers76 REAL(wp) :: zfact, z1_e3t, zdep ! local scalar78 INTEGER :: ji, jj, jk, jn ! dummy loop indices 79 INTEGER :: ikt, ikb ! local integers 80 REAL(wp) :: zfact, z1_e3t, zdep, ztim ! local scalar 77 81 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds 78 82 !!---------------------------------------------------------------------- … … 122 126 DO jj = 2, jpj 123 127 DO ji = fs_2, fs_jpim1 ! vector opt. 124 sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj) ! non solar heat flux 128 IF ( ll_wd ) THEN ! If near WAD point limit the flux for now 129 IF ( sshn(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 130 sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj) ! non solar heat flux 131 ELSE IF ( sshn(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 132 sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj) * (tanh(5._wp*( ( sshn(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1)) ) 133 ELSE 134 sbc_tsc(ji,jj,jp_tem) = 0._wp 135 ENDIF 136 ELSE 137 sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj) ! non solar heat flux 138 ENDIF 139 125 140 sbc_tsc(ji,jj,jp_sal) = r1_rau0 * sfx(ji,jj) ! salt flux due to freezing/melting 126 141 END DO … … 208 223 IF( iom_use('rnf_x_sss') ) CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) ) ! runoff term on sss 209 224 225 #if defined key_asminc 226 ! 227 !---------------------------------------- 228 ! Assmilation effects 229 !---------------------------------------- 230 ! 231 IF( ln_sshinc ) THEN ! input of heat and salt due to assimilation 232 ! 233 IF( ln_linssh ) THEN 234 DO jj = 2, jpj 235 DO ji = fs_2, fs_jpim1 236 ztim = ssh_iau(ji,jj) / e3t_n(ji,jj,1) 237 tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + tsn(ji,jj,1,jp_tem) * ztim 238 tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + tsn(ji,jj,1,jp_sal) * ztim 239 END DO 240 END DO 241 ELSE 242 DO jj = 2, jpj 243 DO ji = fs_2, fs_jpim1 244 ztim = ssh_iau(ji,jj) / ( ht_n(ji,jj) + 1. - ssmask(ji, jj) ) 245 tsa(ji,jj,:,jp_tem) = tsa(ji,jj,:,jp_tem) + tsn(ji,jj,:,jp_tem) * ztim 246 tsa(ji,jj,:,jp_sal) = tsa(ji,jj,:,jp_sal) + tsn(ji,jj,:,jp_sal) * ztim 247 END DO 248 END DO 249 ENDIF 250 ! 251 ENDIF 252 ! 253 #endif 254 210 255 ! 211 256 !---------------------------------------- -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/stpctl.F90
r7852 r8992 22 22 USE lib_mpp ! distributed memory computing 23 23 USE lib_fortran ! Fortran routines library 24 USE wet_dry, ONLY: ll_wd, ssh_ref ! reference depth for negative bathy 24 25 25 26 IMPLICIT NONE … … 150 151 DO jj = 1, jpj 151 152 DO ji = 1, jpi 152 IF( tmask(ji,jj,1) == 1) zsshmax = MAX( zsshmax, ABS(sshn(ji,jj)) ) 153 IF( ll_wd ) THEN 154 IF( tmask(ji,jj,1) == 1) zsshmax = MAX( zsshmax, ABS(sshn(ji,jj)+ssh_ref) ) 155 ELSE 156 IF( tmask(ji,jj,1) == 1) zsshmax = MAX( zsshmax, ABS(sshn(ji,jj)) ) 157 ENDIF 153 158 END DO 154 159 END DO
Note: See TracChangeset
for help on using the changeset viewer.