!! Fortran version of different diagnostics ! L. Fita. LMD May 2016 ! gfortran module_generic.o module_ForDiagnosticsVars.o -c module_ForDiagnostics.F90 ! ! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90 MODULE module_ForDiagnostics USE module_definitions USE module_generic USE module_ForDiagnosticsVars CONTAINS !!!!!!! Calculations ! compute_cape_afwa4D: Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute ! CAPE, CIN, ZLFC, PLFC, LI ! compute_cellbnds: Subroutine to compute cellboundaries using wind-staggered lon, lats as ! intersection of their related parallels and meridians ! compute_cellbndsreg: Subroutine to compute cellboundaries using lon, lat from a reglar lon/lat ! projection as intersection of their related parallels and meridians ! compute_cllmh4D3: Computation of low, medium and high cloudiness from a 4D CLDFRA and pressure being ! 3rd dimension the z-dim ! compute_cllmh3D3: Computation of low, medium and high cloudiness from a 3D CLDFRA and pressure being ! 3rd dimension the z-dim ! compute_cllmh: Computation of low, medium and high cloudiness ! compute_clt4D3: Computation of total cloudiness from a 4D CLDFRA being 3rd dimension the z-dim ! compute_clt3D3: Computation of total cloudiness from a 3D CLDFRA being 3rd dimension the z-dim ! compute_clt: Computation of total cloudiness ! compute_fog_K84: Computation of fog and visibility following Kunkel, (1984) ! compute_fog_RUC: Computation of fog and visibility following RUC method Smirnova, (2000) ! compute_fog_FRAML50: fog and visibility following Gultepe and Milbrandt, (2010) ! compute_massvertint1D: Subroutine to vertically integrate a 1D variable in eta vertical coordinates ! compute_psl_ecmwf: Compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa] ! compute_range_faces: Subroutine to compute faces [uphill, valleys, downhill] of a mountain range along a given face ! compute_vertint1D: Subroutine to vertically integrate a 1D variable in any vertical coordinates ! compute_zint4D: Subroutine to vertically integrate a 4D variable in any vertical coordinates ! compute_zmla_generic4D: Subroutine to compute pbl-height following a generic method ! compute_zwind4D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology ! compute_zwind_log4D: Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology ! compute_zwindMCO3D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodolog !!! ! Calculations !!! SUBROUTINE compute_cllmh4D2(cldfra4D, pres4D, cllmh4D2, d1, d2, d3, d4) ! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA and pressure ! where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> cllmh(3,d1,d3,d4) 1: low, 2: medium, 3: high ! It should be properly done via an 'INTERFACE', but... IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4 REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: cldfra4D, pres4D REAL(r_k), DIMENSION(3,d1,d3,d4), INTENT(out) :: cllmh4D2 ! Local INTEGER :: i,j,k, zdim, Ndim !!!!!!! Variables ! cldfra4D: 4D cloud fraction values [1] ! pres4D: 4D pressure values [Pa] ! Ndim: number of dimensions of the input data ! d[1-4]: dimensions of 'cldfra4D' ! zdim: number of the vertical-dimension within the matrix ! cltlmh4D2: low, medium, high cloudiness for the 4D cldfra and d2 being 'zdim' fname = 'compute_cllmh4D2' zdim = 2 Ndim = 4 DO i=1, d1 DO j=1, d3 DO k=1, d4 cllmh4D2(:,i,j,k) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2) END DO END DO END DO RETURN END SUBROUTINE compute_cllmh4D2 SUBROUTINE compute_cllmh3D1(cldfra3D, pres3D, cllmh3D1, d1, d2, d3) ! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA and pressure ! where zdim is the 1st dimension (thus, cldfra3D(d1,d2,d3) --> cllmh(3,d2,d3) 1: low, 2: medium, 3: high ! It should be properly done via an 'INTERFACE', but... IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3 REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: cldfra3D, pres3D REAL(r_k), DIMENSION(3,d2,d3), INTENT(out) :: cllmh3D1 ! Local INTEGER :: i,j,k, zdim, Ndim !!!!!!! Variables ! cldfra3D: 3D cloud fraction values [1] ! pres3D: 3D pressure values [Pa] ! Ndim: number of dimensions of the input data ! d[1-3]: dimensions of 'cldfra3D' ! zdim: number of the vertical-dimension within the matrix ! cltlmh3D1: low, medium, high cloudiness for the 3D cldfra and d1 being 'zdim' fname = 'compute_cllmh3D1' zdim = 1 Ndim = 3 DO i=1, d1 DO j=1, d2 cllmh3D1(:,i,j) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1) END DO END DO RETURN END SUBROUTINE compute_cllmh3D1 SUBROUTINE compute_cllmh(cldfra1D, cldfra2D, cldfra3D, cldfra4D, pres1D, pres2D, pres3D, pres4D, & Ndim, zdim, cllmh1D, cllmh2D1, cllmh2D2, cllmh3D1, cllmh3D2, cllmh3D3, cllmh4D1, cllmh4D2, & cllmh4D3, cllmh4D4, d1, d2, d3, d4) ! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ IMPLICIT NONE INTEGER, INTENT(in) :: Ndim, d1, d2, d3, d4, zdim REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in) :: cldfra1D, pres1D REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in) :: cldfra2D, pres2D REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D, pres3D REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL, & INTENT(in) :: cldfra4D, pres4D REAL(r_k), DIMENSION(3), OPTIONAL, INTENT(out) :: cllmh1D REAL(r_k), DIMENSION(d1,3), OPTIONAL, INTENT(out) :: cllmh2D1 REAL(r_k), DIMENSION(d2,3), OPTIONAL, INTENT(out) :: cllmh2D2 REAL(r_k), DIMENSION(d2,d3,3), OPTIONAL, INTENT(out) :: cllmh3D1 REAL(r_k), DIMENSION(d1,d3,3), OPTIONAL, INTENT(out) :: cllmh3D2 REAL(r_k), DIMENSION(d1,d2,3), OPTIONAL, INTENT(out) :: cllmh3D3 REAL(r_k), DIMENSION(d2,d3,d4,3), OPTIONAL, & INTENT(out) :: cllmh4D1 REAL(r_k), DIMENSION(d1,d3,d4,3), OPTIONAL, & INTENT(out) :: cllmh4D2 REAL(r_k), DIMENSION(d1,d2,d4,3), OPTIONAL, & INTENT(out) :: cllmh4D3 REAL(r_k), DIMENSION(d1,d2,d3,3), OPTIONAL, & INTENT(out) :: cllmh4D4 ! Local INTEGER :: i,j,k !!!!!!! Variables ! cldfra[1-4]D: cloud fraction values [1] ! pres[1-4]D: pressure values [Pa] ! Ndim: number of dimensions of the input data ! d[1-4]: dimensions of 'cldfra' ! zdim: number of the vertical-dimension within the matrix ! cllmh1D: low, medium and high cloudiness for the 1D cldfra ! cllmh2D1: low, medium and high cloudiness for the 2D cldfra and d1 being 'zdim' ! cllmh2D2: low, medium and high cloudiness for the 2D cldfra and d2 being 'zdim' ! cllmh3D1: low, medium and high cloudiness for the 3D cldfra and d1 being 'zdim' ! cllmh3D2: low, medium and high cloudiness for the 3D cldfra and d2 being 'zdim' ! cllmh3D3: low, medium and high cloudiness for the 3D cldfra and d3 being 'zdim' ! cllmh4D1: low, medium and high cloudiness for the 4D cldfra and d1 being 'zdim' ! cllmh4D2: low, medium and high cloudiness for the 4D cldfra and d2 being 'zdim' ! cllmh4D3: low, medium and high cloudiness for the 4D cldfra and d3 being 'zdim' ! cllmh4D4: low, medium and high cloudiness for the 4D cldfra and d4 being 'zdim' fname = 'compute_cllmh' SELECT CASE (Ndim) CASE (1) cllmh1D = var_cllmh(cldfra1D, pres1D, d1) CASE (2) IF (zdim == 1) THEN DO i=1, d2 cllmh2D1(i,:) = var_cllmh(cldfra2D(:,i), pres2D(:,i), d1) END DO ELSE IF (zdim == 2) THEN DO i=1, d1 cllmh2D2(i,:) = var_cllmh(cldfra2D(:,i), pres2D(i,:), d2) END DO ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!' PRINT *,' accepted values: 1,2' STOP END IF CASE (3) IF (zdim == 1) THEN DO i=1, d2 DO j=1, d3 cllmh3D1(i,j,:) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1) END DO END DO ELSE IF (zdim == 2) THEN DO i=1, d1 DO j=1, d3 cllmh3D2(i,j,:) = var_cllmh(cldfra3D(i,:,j), pres3D(i,:,j), d2) END DO END DO ELSE IF (zdim == 3) THEN DO i=1, d1 DO j=1, d2 cllmh3D3(i,j,:) = var_cllmh(cldfra3D(i,j,:), pres3D(i,j,:), d3) END DO END DO ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!' PRINT *,' accepted values: 1,2,3' STOP END IF CASE (4) IF (zdim == 1) THEN DO i=1, d2 DO j=1, d3 DO k=1, d4 cllmh4D1(i,j,k,:) = var_cllmh(cldfra4D(:,i,j,k), pres4D(:,i,j,k), d1) END DO END DO END DO ELSE IF (zdim == 2) THEN DO i=1, d1 DO j=1, d3 DO k=1, d4 cllmh4D2(i,j,k,:) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2) END DO END DO END DO ELSE IF (zdim == 3) THEN DO i=1, d2 DO j=1, d3 DO k=1, d4 cllmh4D3(i,j,k,:) = var_cllmh(cldfra4D(i,j,:,k), pres4D(i,j,:,k), d3) END DO END DO END DO ELSE IF (zdim == 4) THEN DO i=1, d1 DO j=1, d2 DO k=1, d3 cllmh4D4(i,j,k,:) = var_cllmh(cldfra4D(i,j,k,:), pres4D(i,j,k,:), d4) END DO END DO END DO ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!' PRINT *,' accepted values: 1,2,3,4' STOP END IF CASE DEFAULT PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!' STOP END SELECT RETURN END SUBROUTINE compute_cllmh SUBROUTINE compute_clt4D2(cldfra4D, clt4D2, d1, d2, d3, d4) ! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA ! where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> clt(d1,d3,d4) ! It should be properly done via an 'INTERFACE', but... IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4 REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: cldfra4D REAL(r_k), DIMENSION(d1,d3,d4), INTENT(out) :: clt4D2 ! Local INTEGER :: i,j,k, zdim, Ndim !!!!!!! Variables ! cldfra4D: 4D cloud fraction values [1] ! Ndim: number of dimensions of the input data ! d[1-4]: dimensions of 'cldfra4D' ! zdim: number of the vertical-dimension within the matrix ! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim' fname = 'compute_clt4D2' zdim = 2 Ndim = 4 DO i=1, d1 DO j=1, d3 DO k=1, d4 clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2) END DO END DO END DO RETURN END SUBROUTINE compute_clt4D2 SUBROUTINE compute_clt3D1(cldfra3D, clt3D1, d1, d2, d3) ! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA ! where zdim is the 1st dimension (thus, cldfra4D(d1,d2,d3) --> clt(d2,d3) ! It should be properly done via an 'INTERFACE', but... IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3 REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: cldfra3D REAL(r_k), DIMENSION(d2,d3), INTENT(out) :: clt3D1 ! Local INTEGER :: i,j,k, zdim, Ndim !!!!!!! Variables ! cldfra3D: 3D cloud fraction values [1] ! Ndim: number of dimensions of the input data ! d[1-3]: dimensions of 'cldfra3D' ! zdim: number of the vertical-dimension within the matrix ! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim' fname = 'compute_clt3D1' zdim = 1 Ndim = 3 DO i=1, d2 DO j=1, d3 clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1) END DO END DO RETURN END SUBROUTINE compute_clt3D1 SUBROUTINE compute_clt(cldfra1D, cldfra2D, cldfra3D, cldfra4D, Ndim, zdim, clt1D, clt2D1, clt2D2, & clt3D1, clt3D2, clt3D3, clt4D1, clt4D2, clt4D3, clt4D4, d1, d2, d3, d4) ! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ IMPLICIT NONE INTEGER, INTENT(in) :: Ndim, d1, d2, d3, d4, zdim REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in) :: cldfra1D REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in) :: cldfra2D REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL, & INTENT(in) :: cldfra4D REAL(r_k), OPTIONAL, INTENT(out) :: clt1D REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(out) :: clt2D1 REAL(r_k), DIMENSION(d2), OPTIONAL, INTENT(out) :: clt2D2 REAL(r_k), DIMENSION(d2,d3), OPTIONAL, INTENT(out) :: clt3D1 REAL(r_k), DIMENSION(d1,d3), OPTIONAL, INTENT(out) :: clt3D2 REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(out) :: clt3D3 REAL(r_k), DIMENSION(d2,d3,d4), OPTIONAL,INTENT(out) :: clt4D1 REAL(r_k), DIMENSION(d1,d3,d4), OPTIONAL,INTENT(out) :: clt4D2 REAL(r_k), DIMENSION(d1,d2,d4), OPTIONAL,INTENT(out) :: clt4D3 REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL,INTENT(out) :: clt4D4 ! Local INTEGER :: i,j,k !!!!!!! Variables ! cldfra[1-4]D: cloud fraction values [1] ! Ndim: number of dimensions of the input data ! d[1-4]: dimensions of 'cldfra' ! zdim: number of the vertical-dimension within the matrix ! clt1D: total cloudiness for the 1D cldfra ! clt2D1: total cloudiness for the 2D cldfra and d1 being 'zdim' ! clt2D2: total cloudiness for the 2D cldfra and d2 being 'zdim' ! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim' ! clt3D2: total cloudiness for the 3D cldfra and d2 being 'zdim' ! clt3D3: total cloudiness for the 3D cldfra and d3 being 'zdim' ! clt4D1: total cloudiness for the 4D cldfra and d1 being 'zdim' ! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim' ! clt4D3: total cloudiness for the 4D cldfra and d3 being 'zdim' ! clt4D4: total cloudiness for the 4D cldfra and d4 being 'zdim' fname = 'compute_clt' SELECT CASE (Ndim) CASE (1) clt1D = var_clt(cldfra1D, d1) CASE (2) IF (zdim == 1) THEN DO i=1, d2 clt2D1(i) = var_clt(cldfra2D(:,i), d1) END DO ELSE IF (zdim == 2) THEN DO i=1, d1 clt2D2(i) = var_clt(cldfra2D(:,i), d2) END DO ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!' PRINT *,' accepted values: 1,2' STOP END IF CASE (3) IF (zdim == 1) THEN DO i=1, d2 DO j=1, d3 clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1) END DO END DO ELSE IF (zdim == 2) THEN DO i=1, d1 DO j=1, d3 clt3D2(i,j) = var_clt(cldfra3D(i,:,j), d2) END DO END DO ELSE IF (zdim == 3) THEN DO i=1, d1 DO j=1, d2 clt3D3(i,j) = var_clt(cldfra3D(i,j,:), d3) END DO END DO ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!' PRINT *,' accepted values: 1,2,3' STOP END IF CASE (4) IF (zdim == 1) THEN DO i=1, d2 DO j=1, d3 DO k=1, d4 clt4D1(i,j,k) = var_clt(cldfra4D(:,i,j,k), d1) END DO END DO END DO ELSE IF (zdim == 2) THEN DO i=1, d1 DO j=1, d3 DO k=1, d4 clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2) END DO END DO END DO ELSE IF (zdim == 3) THEN DO i=1, d2 DO j=1, d3 DO k=1, d4 clt4D3(i,j,k) = var_clt(cldfra4D(i,j,:,k), d3) END DO END DO END DO ELSE IF (zdim == 4) THEN DO i=1, d1 DO j=1, d2 DO k=1, d3 clt4D4(i,j,k) = var_clt(cldfra4D(i,j,k,:), d4) END DO END DO END DO ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!' PRINT *,' accepted values: 1,2,3,4' STOP END IF CASE DEFAULT PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!' STOP END SELECT RETURN END SUBROUTINE compute_clt SUBROUTINE compute_massvertint1D(var, mutot, dz, deta, integral) ! Subroutine to vertically integrate a 1D variable in eta vertical coordinates IMPLICIT NONE INTEGER, INTENT(in) :: dz REAL(r_k), INTENT(in) :: mutot REAL(r_k), DIMENSION(dz), INTENT(in) :: var, deta REAL(r_k), INTENT(out) :: integral ! Local INTEGER :: k !!!!!!! Variables ! var: vertical variable to integrate (assuming kgkg-1) ! mutot: total dry-air mass in column ! dz: vertical dimension ! deta: eta-levels difference between full eta-layers fname = 'compute_massvertint1D' ! integral=0. ! DO k=1,dz ! integral = integral + var(k)*deta(k) ! END DO integral = SUM(var*deta) integral=integral*mutot/g RETURN END SUBROUTINE compute_massvertint1D SUBROUTINE compute_zint4D(var4D, dlev, zweight, d1, d2, d3, d4, int3D) ! Subroutine to vertically integrate a 4D variable in any vertical coordinates IMPLICIT NONE INTEGER, INTENT(in) :: d1,d2,d3,d4 REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: var4D, dlev, zweight REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: int3D ! Local INTEGER :: i,j,l !!!!!!! Variables ! var4D: vertical variable to integrate ! dlev: height of layers ! zweight: weight for each level to be applied (=1. for no effect) fname = 'compute_zint4D' DO i=1,d1 DO j=1,d2 DO l=1,d4 CALL compute_vertint1D(var4D(i,j,:,l),d3, dlev(i,j,:,l), zweight(i,j,:,l), & int3D(i,j,l)) END DO END DO END DO RETURN END SUBROUTINE compute_zint4D SUBROUTINE compute_vertint1D(var, dz, deta, zweight, integral) ! Subroutine to vertically integrate a 1D variable in any vertical coordinates IMPLICIT NONE INTEGER, INTENT(in) :: dz REAL(r_k), DIMENSION(dz), INTENT(in) :: var, deta, zweight REAL(r_k), INTENT(out) :: integral ! Local INTEGER :: k !!!!!!! Variables ! var: vertical variable to integrate ! dz: vertical dimension ! deta: eta-levels difference between layers ! zweight: weight for each level to be applied (=1. for no effect) fname = 'compute_vertint1D' ! integral=0. ! DO k=1,dz ! integral = integral + var(k)*deta(k) ! END DO integral = SUM(var*deta*zweight) RETURN END SUBROUTINE compute_vertint1D SUBROUTINE compute_cape_afwa4D(ta, hur, press, zg, hgt, cape, cin, zlfc, plfc, li, parcelmethod, & d1, d2, d3, d4) ! Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute CAPE, CIN, ZLFC, PLFC, LI IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4, parcelmethod REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ta, hur, press, zg REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: cape, cin, zlfc, plfc, li ! Local INTEGER :: i, j, it INTEGER :: ofunc !!!!!!! Variables ! ta: air temperature [K] ! hur: relative humidity [%] ! press: air pressure [Pa] ! zg: geopotential height [gpm] ! hgt: topographical height [m] ! cape: Convective available potential energy [Jkg-1] ! cin: Convective inhibition [Jkg-1] ! zlfc: height at the Level of free convection [m] ! plfc: pressure at the Level of free convection [Pa] ! li: lifted index [1] ! parcelmethod: ! Most Unstable = 1 (default) ! Mean layer = 2 ! Surface based = 3 fname = 'compute_cape_afwa4D' DO i=1, d1 DO j=1, d2 DO it=1, d4 ofunc = var_cape_afwa1D(d3, ta(i,j,:,it), hur(i,j,:,it), press(i,j,:,it), zg(i,j,:,it), & 1, cape(i,j,it), cin(i,j,it), zlfc(i,j,it), plfc(i,j,it), li(i,j,it), parcelmethod) IF (zlfc(i,j,it) /= -1.) zlfc(i,j,it) = zlfc(i,j,it) - hgt(i,j) END DO END DO END DO RETURN END SUBROUTINE compute_cape_afwa4D SUBROUTINE compute_psl_ecmwf(ps, hgt, T, press, unpress, psl, d1, d2, d4) ! Subroutine to compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa] IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d4 REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: ps, T, press, unpress REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: psl ! Local INTEGER :: i, j, it !!!!!!! Variables ! ps: surface pressure [Pa] ! hgt: terrain height [m] ! T: temperature at first half-mass level [K] ! press: pressure at first full levels [Pa] ! unpress: pressure at first mass (half) levels [Pa] ! psl: sea-level pressure [Pa] fname = 'compute_psl_ecmwf' DO i=1, d1 DO j=1, d2 DO it=1, d4 CALL var_psl_ecmwf(ps(i,j,it), hgt(i,j), T(i,j,it), unpress(i,j,it), press(i,j,it), & psl(i,j,it)) END DO END DO END DO RETURN END SUBROUTINE compute_psl_ecmwf SUBROUTINE compute_zmla_generic4D(tpot, qratio, z, hgt, zmla3D, d1, d2, d3, d4) ! Subroutine to compute pbl-height following a generic method ! from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim. ! applied also in Garcia-Diez et al., 2013, QJRMS ! where ! "The technique identifies the ML height as a threshold increase of potential temperature from ! its minimum value within the boundary layer." ! here applied similarly to Garcia-Diez et al. where ! zmla = "...first level where potential temperature exceeds the minimum potential temperature ! reached in the mixed layer by more than 1.5 K" IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4 REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: tpot, qratio, z REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: zmla3D ! Local INTEGER :: i, j, it !!!!!!! Variables ! tpot: potential air temperature [K] ! qratio: water vapour mixing ratio [kgkg-1] ! z: height above sea level [m] ! hgt: terrain height [m] ! zmla3D: boundary layer height from surface [m] fname = 'compute_zmla_generic4D' DO i=1, d1 DO j=1, d2 DO it=1, d4 CALL var_zmla_generic(d3, qratio(i,j,:,it), tpot(i,j,:,it), z(i,j,:,it), hgt(i,j), & zmla3D(i,j,it)) END DO END DO END DO RETURN END SUBROUTINE compute_zmla_generic4D SUBROUTINE compute_zwind4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4) ! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4 REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ua, va, z REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: uas, vas REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa REAL(r_k), INTENT(in) :: zextrap REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: uaz, vaz ! Local INTEGER :: i, j, it !!!!!!! Variables ! tpot: potential air temperature [K] ! qratio: water vapour mixing ratio [kgkg-1] ! z: height above surface [m] ! sina, cosa: local sine and cosine of map rotation [1.] ! zmla3D: boundary layer height from surface [m] fname = 'compute_zwind4D' DO i=1, d1 DO j=1, d2 DO it=1, d4 CALL var_zwind(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it), & sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it)) END DO END DO END DO RETURN END SUBROUTINE compute_zwind4D SUBROUTINE compute_zwind_log4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4) ! Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4 REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ua, va, z REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: uas, vas REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa REAL(r_k), INTENT(in) :: zextrap REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: uaz, vaz ! Local INTEGER :: i, j, it !!!!!!! Variables ! tpot: potential air temperature [K] ! qratio: water vapour mixing ratio [kgkg-1] ! z: height above surface [m] ! sina, cosa: local sine and cosine of map rotation [1.] ! zmla3D: boundary layer height from surface [m] fname = 'compute_zwind_log4D' DO i=1, d1 DO j=1, d2 DO it=1, d4 CALL var_zwind_log(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it), & sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it)) END DO END DO END DO RETURN END SUBROUTINE compute_zwind_log4D SUBROUTINE compute_zwindMO3D(d1, d2, d3, ust, znt, rmol, uas, vas, sina, cosa, newz, uznew, vznew) ! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology ! NOTE: only usefull for newz < 80. m IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3 REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: ust, znt, rmol REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: uas, vas REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa REAL(r_k), INTENT(in) :: newz REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: uznew, vznew ! Local INTEGER :: i, j, it !!!!!!! Variables ! ust: u* in similarity theory [ms-1] ! znt: thermal time-varying roughness length [m] ! rmol: Inverse of the Obukhov length [m-1] ! uas: x-component 10-m wind speed [ms-1] ! vas: y-component 10-m wind speed [ms-1] ! sina, cosa: local sine and cosine of map rotation [1.] fname = 'compute_zwindMO3D' DO i=1, d1 DO j=1, d2 DO it=1, d3 CALL var_zwind_MOtheor(ust(i,j,it), znt(i,j,it), rmol(i,j,it), uas(i,j,it), vas(i,j,it), & sina(i,j), cosa(i,j), newz, uznew(i,j,it), vznew(i,j,it)) END DO END DO END DO RETURN END SUBROUTINE compute_zwindMO3D SUBROUTINE compute_potevap_orPM3D(d1, d2, d3, rho1, ust, uas, vas, tas, ps, qv1, potevap) ! Subroutine to compute potential evapotranspiration Penman-Monteith formulation implemented in ! ORCHIDEE in src_sechiba/enerbil.f90 IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3 REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: rho1, ust, uas, vas, tas, ps, qv1 REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: potevap ! Local INTEGER :: i, j, it !!!!!!! Variables ! rho1: atsmophere density at the first layer [kgm-3] ! ust: u* in similarity theory [ms-1] ! uas: x-component 10-m wind speed [ms-1] ! vas: y-component 10-m wind speed [ms-1] ! tas: 2-m atmosphere temperature [K] ! ps: surface pressure [Pa] ! qv1: 1st layer atmospheric mixing ratio [kgkg-1] ! potevap: potential evapo transpiration [kgm-2s-1] fname = 'compute_potevap_orPM3D' DO i=1, d1 DO j=1, d2 DO it=1, d3 CALL var_potevap_orPM(rho1(i,j,it), ust(i,j,it), uas(i,j,it), vas(i,j,it), tas(i,j,it), & ps(i,j,it), qv1(i,j,it), potevap(i,j,it)) END DO END DO END DO RETURN END SUBROUTINE compute_potevap_orPM3D SUBROUTINE compute_fog_K84(d1, d2, d3, qc, qi, fog, vis) ! Subroutine to compute fog: qcloud + qice /= 0. ! And visibility following Kunkel, B. A., (1984): Parameterization of droplet terminal velocity and ! extinction coefficient in fog models. J. Climate Appl. Meteor., 23, 34–41. IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3 REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: qc, qi INTEGER, DIMENSION(d1,d2,d3), INTENT(out) :: fog REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: vis ! Local INTEGER :: i, j, it !!!!!!! Variables ! qc: cloud mixing ratio [kgkg-1] ! qi, ice mixing ratio [kgkg-1] ! fog: presence of fog (1: yes, 0: no) ! vis: visibility within fog [km] fname = 'compute_fog_K84' DO i=1, d1 DO j=1, d2 DO it=1, d3 CALL var_fog_K84(qc(i,j,it), qi(i,j,it), fog(i,j,it), vis(i,j,it)) END DO END DO END DO RETURN END SUBROUTINE compute_fog_K84 SUBROUTINE compute_fog_RUC(d1, d2, d3, qv, ta, pres, fog, vis) ! Subroutine to compute fog: qcloud + qice /= 0. ! And visibility following RUC method Smirnova, T. G., S. G. Benjamin, and J. M. Brown, 2000: Case ! study verification of RUC/MAPS fog and visibility forecasts. Preprints, 9 th Conference on ! Aviation, Range, and Aerospace Meteorlogy, AMS, Orlando, FL, Sep. 2000. Paper#2.3, 6 pp. IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3 REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: qv, ta, pres INTEGER, DIMENSION(d1,d2,d3), INTENT(out) :: fog REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: vis ! Local INTEGER :: i, j, it !!!!!!! Variables ! qv: water vapor mixing ratio [kgkg-1] ! ta: temperature [K] ! pres: pressure [Pa] ! fog: presence of fog (1: yes, 0: no) ! vis: visibility within fog [km] fname = 'compute_fog_RUC' DO i=1, d1 DO j=1, d2 DO it=1, d3 CALL var_fog_RUC(qv(i,j,it), ta(i,j,it), pres(i,j,it), fog(i,j,it), vis(i,j,it)) END DO END DO END DO RETURN END SUBROUTINE compute_fog_RUC SUBROUTINE compute_fog_FRAML50(d1, d2, d3, qv, ta, pres, fog, vis) ! Subroutine to compute fog (vis < 1 km) and visibility following ! Gultepe, I. and J.A. Milbrandt, 2010: Probabilistic Parameterizations of Visibility Using ! Observations of Rain Precipitation Rate, Relative Humidity, and Visibility. J. Appl. Meteor. ! Climatol., 49, 36-46, https://doi.org/10.1175/2009JAMC1927.1 ! Interest is focused on a 'general' fog/visibilty approach, thus the fit at 50 % of probability is ! chosen ! Effects from precipitation are not considered IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3 REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: qv, ta, pres INTEGER, DIMENSION(d1,d2,d3), INTENT(out) :: fog REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: vis ! Local INTEGER :: i, j, it !!!!!!! Variables ! qv: mixing ratio in [kgkg-1] ! ta: temperature [K] ! pres: pressure field [Pa] ! fog: presence of fog (1: yes, 0: no) ! vis: visibility within fog [km] fname = 'compute_fog_FRAML50' DO i=1, d1 DO j=1, d2 DO it=1, d3 CALL var_fog_FRAML50(qv(i,j,it), ta(i,j,it), pres(i,j,it), fog(i,j,it), vis(i,j,it)) END DO END DO END DO RETURN END SUBROUTINE compute_fog_FRAML50 SUBROUTINE compute_range_faces(d1, d2, lon, lat, hgt, xdist, ydist, dist, face, dsfilt, dsnewrange, & hvalrng, hgtmax, pthgtmax, derivhgt, peaks, valleys, origfaces, filtfaces, ranges, rangeshgtmax, & ptrangeshgtmax) ! Subroutine to compute faces [uphill, valleys, downhill] of a mountain range along a given face IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2 REAL(r_k), INTENT(in) :: dsfilt, dsnewrange, hvalrng REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: lon, lat, hgt, xdist, ydist, dist CHARACTER(len=*) :: face REAL(r_k), DIMENSION(d1,d2), INTENT(out) :: derivhgt, hgtmax, rangeshgtmax INTEGER, DIMENSION(d1,d2), INTENT(out) :: pthgtmax, origfaces, filtfaces, peaks, & valleys, ranges, ptrangeshgtmax ! Local INTEGER :: i, j INTEGER :: pthgtmax1, Npeaks, Nvalleys, Nranges REAL(r_k) :: hgtmax1 INTEGER, DIMENSION(d1) :: ipeaks1, ivalleys1, irangeshgtmax1 INTEGER, DIMENSION(d2) :: ipeaks2, ivalleys2, irangeshgtmax2 REAL(r_k), DIMENSION(d1) :: rangeshgtmax1 REAL(r_k), DIMENSION(d2) :: rangeshgtmax2 INTEGER, DIMENSION(2,d1) :: ranges1 INTEGER, DIMENSION(2,d2) :: ranges2 !!!!!!! Variables ! lon: longitude [degrees east] ! lat: latitude [degrees north] ! hgt: topograpical height [m] ! face: which face (axis along which produce slices) to use to compute the faces: WE, SN ! dsfilt: distance to filter orography smaller scale of it [m] ! dsnewrange: distance to start a new mountain range [m] ! hvalrng: maximum height of a valley to mark change of range [m] ! hgtmax: maximum height of the face [m] ! pthgtmax: grid point of the maximum height [1] ! derivhgt: topograpic derivative along axis [m deg-1] ! peaks: peak point ! valleys: valley point ! origfaces: original faces (-1, downhill; 0: valley; 1: uphill) ! filtfaces: filtered faces (-1, downhill; 0: valley; 1: uphill) ! ranges: number of range ! rangeshgtmax: maximum height for each individual range [m] ! ptrangeshgtmax: grid point maximum height for each individual range [1] fname = 'compute_range_faces' peaks = 0 valleys = 0 pthgtmax = 0 rangeshgtmax = fillVal64 IF (TRIM(face) == 'WE') THEN DO j=1, d2 !PRINT *,'Lluis:', j-1, '***' CALL var_range_faces(d1, lon(:,j), lat(:,j), hgt(:,j), xdist(:,j), dsfilt, & dsnewrange, hvalrng, hgtmax1, pthgtmax1, derivhgt(:,j), Npeaks, ipeaks1, & Nvalleys, ivalleys1, origfaces(:,j), filtfaces(:,j), Nranges, ranges1, & rangeshgtmax1, irangeshgtmax1) hgtmax(:,j) = hgtmax1 pthgtmax(pthgtmax1,j) = 1 DO i=1, Npeaks peaks(ipeaks1(i),j) = 1 END DO DO i=1, Nvalleys valleys(ivalleys1(i),j) = 1 END DO DO i=1, Nranges ranges(ranges1(1,i):ranges1(2,i),j) = i rangeshgtmax(ranges1(1,i):ranges1(2,i),j) = rangeshgtmax1(i) ptrangeshgtmax(irangeshgtmax1(i),j) = 1 END DO END DO ELSE IF (TRIM(face) == 'SN') THEN DO i=1, d1 CALL var_range_faces(d2, lon(i,:), lat(i,:), hgt(i,:), ydist(i,:), dsfilt, & dsnewrange, hvalrng, hgtmax1, pthgtmax1, derivhgt(i,:), Npeaks, ipeaks2, & Nvalleys, ivalleys2, origfaces(i,:), filtfaces(i,:), Nranges, ranges2, & rangeshgtmax2, irangeshgtmax2) hgtmax(i,:) = hgtmax1 pthgtmax(i,pthgtmax1) = 1 DO j=1, Npeaks peaks(i,ipeaks2(j)) = 1 END DO DO j=1, Nvalleys valleys(i,ivalleys2(j)) = 1 END DO DO j=1, Nranges ranges(i,ranges2(1,j):ranges2(2,j)) = j rangeshgtmax(i,ranges2(1,j):ranges2(2,j)) = rangeshgtmax2(j) ptrangeshgtmax(i,irangeshgtmax2(j)) = 1 END DO END DO ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ": wrong face: '" // TRIM(face) // "' !!" PRINT *,' accepted ones: WE, SN' STOP END IF RETURN END SUBROUTINE compute_range_faces SUBROUTINE compute_cellbnds(dx, dy, sdx, sdy, ulon, ulat, vlon, vlat, xbnds, ybnds) ! Subroutine to compute cellboundaries using wind-staggered lon, lats as intersection of their related ! parallels and meridians IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, sdx, sdy REAL(r_k), DIMENSION(sdx, dy), INTENT(in) :: ulon, ulat REAL(r_k), DIMENSION(dx, sdy), INTENT(in) :: vlon, vlat REAL(r_k), DIMENSION(dx, dy, 4), INTENT(out) :: xbnds, ybnds ! Local INTEGER :: i,j,iv INTEGER :: ix,ex,iy,ey REAL(r_k) :: tmpval1, tmpval2 CHARACTER(len=2), DIMENSION(4) :: Svertex INTEGER, DIMENSION(4,2,2,2) :: indices REAL(r_k), DIMENSION(2) :: ptintsct REAL(r_k), DIMENSION(2,2) :: merid, paral LOGICAL :: intsct !!!!!!! Variables ! dx, dy: un-staggered dimensions ! sdx, sdy: staggered dimensions ! ulon, ulat: x-wind staggered longitudes and latitudes ! vlon, vlat: y-wind staggered longitudes and latitudes ! xbnds, ybnds: x and y cell boundaries fname = 'compute_cellbnds' ! Indices to use indices[SW/NW/NE/SE, m/p, x/y, i/e] Svertex = (/ 'SW', 'NW', 'NE', 'SE' /) ! SW indices(1,1,1,1) = 0 indices(1,1,1,2) = 0 indices(1,1,2,1) = -1 indices(1,1,2,2) = 0 indices(1,2,1,1) = -1 indices(1,2,1,2) = 0 indices(1,2,2,1) = -1 indices(1,2,2,2) = -1 ! NW indices(2,1,1,1) = 0 indices(2,1,1,2) = 0 indices(2,1,2,1) = 0 indices(2,1,2,2) = 1 indices(2,2,1,1) = -1 indices(2,2,1,2) = 0 indices(2,2,2,1) = 1 indices(2,2,2,2) = 1 ! NE indices(3,1,1,1) = 1 indices(3,1,1,2) = 1 indices(3,1,2,1) = 0 indices(3,1,2,2) = 1 indices(3,2,1,1) = 0 indices(3,2,1,2) = 1 indices(3,2,2,1) = 1 indices(3,2,2,2) = 1 ! SE indices(4,1,1,1) = 1 indices(4,1,1,2) = 1 indices(4,1,2,1) = -1 indices(4,1,2,2) = 0 indices(4,2,1,1) = 0 indices(4,2,1,2) = 1 indices(4,2,2,1) = -1 indices(4,2,2,2) = -1 DO i=2,dx-1 DO j=1,dy DO iv=1,4 ix = MAX(i+indices(iv,1,1,1),1) !ex = MIN(i+indices(iv,1,1,2),dx) ex = MAX(i+indices(iv,1,1,2),1) iy = MAX(j+indices(iv,1,2,1),1) ey = MIN(j+indices(iv,1,2,2),dy) merid(1,1) = ulon(ix,iy) merid(1,2) = ulat(ix,iy) merid(2,1) = ulon(ex,ey) merid(2,2) = ulat(ex,ey) ix = MAX(i+indices(iv,2,1,1),1) ex = MIN(i+indices(iv,2,1,2),dx) iy = MAX(j+indices(iv,2,2,1),1) !ey = MIN(i+indices(iv,2,2,2),dy) ey = MAX(j+indices(iv,2,2,2),1) paral(1,1) = vlon(ix,iy) paral(1,2) = vlat(ix,iy) paral(2,1) = vlon(ex,ey) paral(2,2) = vlat(ex,ey) CALL intersection_2Dlines(merid, paral, intsct, ptintsct) IF (.NOT.intsct) THEN msg = 'not intersection found for ' // Svertex(iv) // ' vertex' CALL ErrMsg(msg, fname, -1) END IF xbnds(i,j,iv) = ptintsct(1) ybnds(i,j,iv) = ptintsct(2) END DO END DO END DO ! Dealing with the boundary values i = 1 DO j=1,dy DO iv=1,4 ix = MAX(i+indices(iv,1,1,1),1) !ex = MIN(i+indices(iv,1,1,2),dx) ex = MAX(i+indices(iv,1,1,2),1) iy = MAX(j+indices(iv,1,2,1),1) ey = MIN(j+indices(iv,1,2,2),dy) merid(1,1) = ulon(ix,iy) merid(1,2) = ulat(ix,iy) merid(2,1) = ulon(ex,ey) merid(2,2) = ulat(ex,ey) ix = MAX(i+indices(iv,2,1,1),1) ex = MIN(i+indices(iv,2,1,2),dx) iy = MAX(j+indices(iv,2,2,1),1) !ey = MIN(i+indices(iv,2,2,2),dy) ey = MAX(j+indices(iv,2,2,2),1) IF (iv == 1 .OR. iv == 2) THEN ! Projecting values using dx from next grid point tmpval1 = vlon(2,iy) paral(2,1) = vlon(ex,ey) tmpval2 = tmpval1 - paral(2,1) paral(1,1) = paral(2,1) - tmpval2 tmpval1 = vlat(2,iy) paral(2,2) = vlat(ex,ey) tmpval2 = tmpval1 - paral(2,2) paral(1,2) = paral(2,2) - tmpval2 ELSE paral(1,1) = vlon(ix,iy) paral(1,2) = vlat(ix,iy) paral(2,1) = vlon(ex,ey) paral(2,2) = vlat(ex,ey) END IF CALL intersection_2Dlines(merid, paral, intsct, ptintsct) IF (.NOT.intsct) THEN msg = 'not intersection found for ' // Svertex(iv) // ' vertex' CALL ErrMsg(msg, fname, -1) END IF xbnds(i,j,iv) = ptintsct(1) ybnds(i,j,iv) = ptintsct(2) END DO END DO i = dx DO j=1,dy DO iv=1,4 ix = MAX(i+indices(iv,1,1,1),1) !ex = MIN(i+indices(iv,1,1,2),dx) ex = MAX(i+indices(iv,1,1,2),1) iy = MAX(j+indices(iv,1,2,1),1) ey = MIN(j+indices(iv,1,2,2),dy) merid(1,1) = ulon(ix,iy) merid(1,2) = ulat(ix,iy) merid(2,1) = ulon(ex,ey) merid(2,2) = ulat(ex,ey) ix = MAX(i+indices(iv,2,1,1),1) ex = MIN(i+indices(iv,2,1,2),dx) iy = MAX(j+indices(iv,2,2,1),1) !ey = MIN(i+indices(iv,2,2,2),dy) ey = MAX(j+indices(iv,2,2,2),1) IF (iv == 3 .OR. iv == 4) THEN ! Projecting values using dx from previous grid point tmpval1 = vlon(dx-1,iy) paral(2,1) = vlon(ex,ey) tmpval2 = tmpval1 - paral(2,1) paral(1,1) = paral(2,1) - tmpval2 tmpval1 = vlat(dx-1,iy) paral(2,2) = vlat(ex,ey) tmpval2 = tmpval1 - paral(2,2) paral(1,2) = paral(2,2) - tmpval2 ELSE paral(1,1) = vlon(ix,iy) paral(1,2) = vlat(ix,iy) paral(2,1) = vlon(ex,ey) paral(2,2) = vlat(ex,ey) END IF CALL intersection_2Dlines(merid, paral, intsct, ptintsct) IF (.NOT.intsct) THEN msg = 'not intersection found for ' // Svertex(iv) // ' vertex' CALL ErrMsg(msg, fname, -1) END IF xbnds(i,j,iv) = ptintsct(1) ybnds(i,j,iv) = ptintsct(2) END DO END DO END SUBROUTINE compute_cellbnds SUBROUTINE compute_cellbndsreg(dx, dy, lon, lat, xbnds, ybnds) ! Subroutine to compute cellboundaries using lon, lat from a reglar lon/lat projection as intersection ! of their related parallels and meridians IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx, dy), INTENT(in) :: lon, lat REAL(r_k), DIMENSION(dx, dy, 4), INTENT(out) :: xbnds, ybnds ! Local INTEGER :: i,j,iv INTEGER :: ix,ex,iy,ey CHARACTER(len=2), DIMENSION(4) :: Svertex INTEGER, DIMENSION(4,2,2) :: indices !!!!!!! Variables ! dx, dy: un-staggered dimensions ! lon, lat: longitudes and latitudes ! xbnds, ybnds: x and y cell boundaries fname = 'compute_cellbndsreg' ! Indices to use indices[SW/NW/NE/SE, m/p, x/y, i/e] Svertex = (/ 'SW', 'NW', 'NE', 'SE' /) ! SW indices(1,1,1) = -1 indices(1,1,2) = 0 indices(1,2,1) = -1 indices(1,2,2) = 0 ! NW indices(2,1,1) = -1 indices(2,1,2) = 0 indices(2,2,1) = 0 indices(2,2,2) = 1 ! NE indices(3,1,1) = 0 indices(3,1,2) = 1 indices(3,2,1) = 0 indices(3,2,2) = 1 ! SE indices(4,1,1) = 0 indices(4,1,2) = 1 indices(4,2,1) = -1 indices(4,2,2) = 0 DO i=1,dx DO j=1,dy DO iv=1,4 ix = MAX(i+indices(iv,1,1),1) ix = MIN(ix,dx) ex = MAX(i+indices(iv,1,2),1) ex = MIN(ex,dx) iy = MAX(j+indices(iv,2,1),1) iy = MIN(iy,dy) ey = MAX(j+indices(iv,2,2),1) ey = MIN(ey,dy) xbnds(i,j,iv) = 0.5*(lon(ix,iy) + lon(ex,ey)) ybnds(i,j,iv) = 0.5*(lat(ix,iy) + lat(ex,ey)) END DO END DO END DO END SUBROUTINE SUBROUTINE compute_Koeppen_Geiger_climates(dx, dy, dt, pr, tas, climatesI, climatesS, climlegend) ! Subroutine to compute the Koeppen-Geiger climates after: ! Kottek et al., Meteorologische Zeitschrift, Vol. 15, No. 3, 259-263 IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, dt REAL(r_k), DIMENSION(dx, dy, dt), INTENT(in) :: pr, tas INTEGER, DIMENSION(dy, dt), INTENT(out) :: climatesI CHARACTER(LEN=3), DIMENSION(dy, dt), INTENT(out) :: climatesS CHARACTER(LEN=1000), INTENT(out) :: climlegend ! Local INTEGER :: i,j END SUBROUTINE compute_Koeppen_Geiger_climates END MODULE module_ForDiagnostics