Changeset 12264 for utils/tools/NESTING/src/agrif_readwrite.f90
- Timestamp:
- 2019-12-17T18:32:39+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
utils/tools/NESTING/src/agrif_readwrite.f90
r12253 r12264 545 545 ln_sco = 0 546 546 ln_isfcav = 0 547 IF( partial_steps ) THEN547 !!$ IF( partial_steps ) THEN 548 548 ln_zps = 1 549 549 ln_zco = 0 550 ELSE551 ln_zps = 0552 ln_zco = 1553 ENDIF550 !!$ ELSE 551 !!$ ln_zps = 0 552 !!$ ln_zco = 1 553 !!$ ENDIF 554 554 555 555 CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, & ! input … … 866 866 zmax = gdepw_1d(N) + e3t_1d(N) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(N-1) ) 867 867 zbathy(:,:) = MIN( zmax , zbathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) 868 WHERE( zbathy(:,:) == 0. ) ; mbathy(:,:) = 0 868 WHERE( zbathy(:,:) == 0. ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 869 869 ELSEWHERE ; mbathy(:,:) = N-1 ! ocean : initialize mbathy to the max ocean level 870 870 END WHERE 871 871 872 IF( partial_steps ) THEN872 !!$ IF( partial_steps ) THEN 873 873 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 874 874 ! find the number of ocean levels such that the last level thickness … … 879 879 WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 880 880 END DO 881 ELSE882 DO jk = 1, N883 WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) >= gdepw_1d(jk) ) mbathy(:,:) = jk-1884 END DO885 ENDIF881 !!$ ELSE 882 !!$ DO jk = 1, N 883 !!$ WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) >= gdepw_1d(jk) ) mbathy(:,:) = jk-1 884 !!$ END DO 885 !!$ ENDIF 886 886 887 887 ! Scale factors and depth at T- and W-points … … 897 897 !! IF ( ln_isfcav == 1 ) CALL zgr_isf 898 898 ! 899 IF( partial_steps ) THEN 900 ! Scale factors and depth at T- and W-points 901 ! IF ( ln_isfcav == 0 ) THEN 902 DO jj = 1, ny 903 DO ji = 1, nx 904 ik = mbathy(ji,jj) 905 IF( ik > 0 ) THEN ! ocean point only 906 ! max ocean level case 907 IF( ik == N-1 ) THEN 908 zdepwp = zbathy(ji,jj) 909 ze3tp = zbathy(ji,jj) - gdepw_1d(ik) 910 ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( ze3tp/e3t_1d(ik) ) ) 911 e3t_0(ji,jj,ik ) = ze3tp 912 e3t_0(ji,jj,ik+1) = ze3tp 913 e3w_0(ji,jj,ik ) = ze3wp 914 e3w_0(ji,jj,ik+1) = ze3tp 915 gdepw_0(ji,jj,ik+1) = zdepwp 916 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 917 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 918 ! 919 ELSE ! standard case 920 IF( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = zbathy(ji,jj) 921 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 922 ENDIF 923 !gm Bug? check the gdepw_1d 924 ! ... on ik 925 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 926 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 927 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 928 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 929 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 930 e3w_0(ji,jj,ik) = 0.5 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * gdepw_1d(ik) ) & 931 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 932 ! ... on ik+1 933 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 934 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 935 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 899 ! Scale factors and depth at T- and W-points 900 ! IF ( ln_isfcav == 0 ) THEN 901 DO jj = 1, ny 902 DO ji = 1, nx 903 ik = mbathy(ji,jj) 904 IF( ik > 0 ) THEN ! ocean point only 905 ! max ocean level case 906 IF( ik == N-1 ) THEN 907 zdepwp = zbathy(ji,jj) 908 ze3tp = zbathy(ji,jj) - gdepw_1d(ik) 909 ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( ze3tp/e3t_1d(ik) ) ) 910 e3t_0(ji,jj,ik ) = ze3tp 911 e3t_0(ji,jj,ik+1) = ze3tp 912 e3w_0(ji,jj,ik ) = ze3wp 913 e3w_0(ji,jj,ik+1) = ze3tp 914 gdepw_0(ji,jj,ik+1) = zdepwp 915 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 916 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 917 ! 918 ELSE ! standard case 919 IF( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = zbathy(ji,jj) 920 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 936 921 ENDIF 922 !gm Bug? check the gdepw_1d 923 ! ... on ik 924 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 925 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 926 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 927 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 928 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 929 e3w_0(ji,jj,ik) = 0.5 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * gdepw_1d(ik) ) & 930 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 931 ! ... on ik+1 932 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 933 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 934 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 937 935 ENDIF 938 END DO936 ENDIF 939 937 END DO 940 !941 ! DO jj = 1, ny942 ! DO ji = 1, nx943 ! ik = mbathy(ji,jj)944 ! IF( ik > 0 ) THEN ! ocean point only945 ! e3tp (ji,jj) = e3t_0(ji,jj,ik)946 ! e3wp (ji,jj) = e3w_0(ji,jj,ik)947 ! ENDIF948 ! END DO949 !END DO950 ! END IF951 ENDIF938 END DO 939 ! 940 ! DO jj = 1, ny 941 ! DO ji = 1, nx 942 ! ik = mbathy(ji,jj) 943 ! IF( ik > 0 ) THEN ! ocean point only 944 ! e3tp (ji,jj) = e3t_0(ji,jj,ik) 945 ! e3wp (ji,jj) = e3w_0(ji,jj,ik) 946 ! ENDIF 947 ! END DO 948 ! END DO 949 ! END IF 952 950 ! 953 951 ! Scale factors and depth at U-, V-, UW and VW-points
Note: See TracChangeset
for help on using the changeset viewer.