! $Id: calfis.f90 5182 2024-09-10 14:25:29Z evignon $ SUBROUTINE calfis(lafin, & jD_cur, jH_cur, & pucov, & pvcov, & pteta, & pq, & pmasse, & pps, & pp, & ppk, & pphis, & pphi, & pducov, & pdvcov, & pdteta, & pdq, & flxw, & pdufi, & pdvfi, & pdhfi, & pdqfi, & pdpsfi) ! Auteur : P. Le Van, F. Hourdin ! ......... USE lmdz_infotrac, ONLY: nqtot, tracers USE control_mod, ONLY: planet_type, nsplit_phys USE callphysiq_mod, ONLY: call_physiq USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi USE comvert_mod, ONLY: preff, presnivs USE lmdz_iniprint, ONLY: lunout, prt_level USE lmdz_ssum_scopy, ONLY: scopy, ssum USE lmdz_comgeom2 USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE !======================================================================= ! 1. rearrangement des tableaux et transformation ! variables dynamiques > variables physiques ! 2. calcul des termes physiques ! 3. retransformation des tendances physiques en tendances dynamiques ! remarques: ! ---------- ! - les vents sont donnes dans la physique par leurs composantes ! naturelles. ! - la variable thermodynamique de la physique est une variable ! intensive : T ! pour la dynamique on prend T * ( preff / p(l) ) **kappa ! - les deux seules variables dependant de la geometrie necessaires ! pour la physique sont la latitude pour le rayonnement et ! l'aire de la maille quand on veut integrer une grandeur ! horizontalement. ! - les points de la physique sont les points scalaires de la ! la dynamique; numerotation: ! 1 pour le pole nord ! (jjm-1)*iim pour l'interieur du domaine ! ngridmx pour le pole sud ! ---> ngridmx=2+(jjm-1)*iim ! Input : ! ------- ! pucov covariant zonal velocity ! pvcov covariant meridional velocity ! pteta potential temperature ! pps surface pressure ! pmasse masse d'air dans chaque maille ! pts surface temperature (K) ! callrad clef d'appel au rayonnement ! Output : ! -------- ! pdufi tendency for the natural zonal velocity (ms-1) ! pdvfi tendency for the natural meridional velocity ! pdhfi tendency for the potential temperature ! pdtsfi tendency for the surface temperature ! pdtrad radiative tendencies \ both input ! pfluxrad radiative fluxes / and output !======================================================================= !----------------------------------------------------------------------- ! 0. Declarations : ! ------------------ INTEGER :: ngridmx PARAMETER(ngridmx = 2 + (jjm - 1) * iim - 1 / jjm) ! Arguments : ! ----------- LOGICAL, INTENT(IN) :: lafin ! .TRUE. for the very last CALL to physics REAL, INTENT(IN) :: jD_cur, jH_cur REAL, INTENT(IN) :: pvcov(iip1, jjm, llm) ! covariant meridional velocity REAL, INTENT(IN) :: pucov(iip1, jjp1, llm) ! covariant zonal velocity REAL, INTENT(IN) :: pteta(iip1, jjp1, llm) ! potential temperature REAL, INTENT(IN) :: pmasse(iip1, jjp1, llm) ! mass in each cell ! not used REAL, INTENT(IN) :: pq(iip1, jjp1, llm, nqtot) ! tracers REAL, INTENT(IN) :: pphis(iip1, jjp1) ! surface geopotential REAL, INTENT(IN) :: pphi(iip1, jjp1, llm) ! geopotential REAL, INTENT(IN) :: pdvcov(iip1, jjm, llm) ! dynamical tendency on vcov REAL, INTENT(IN) :: pducov(iip1, jjp1, llm) ! dynamical tendency on ucov REAL, INTENT(IN) :: pdteta(iip1, jjp1, llm) ! dynamical tendency on teta ! NB: pdteta is used only to compute pcvgt which is in fact not used... REAL, INTENT(IN) :: pdq(iip1, jjp1, llm, nqtot) ! dynamical tendency on tracers ! NB: pdq is only used to compute pcvgq which is in fact not used... REAL, INTENT(IN) :: pps(iip1, jjp1) ! surface pressure (Pa) REAL, INTENT(IN) :: pp(iip1, jjp1, llmp1) ! pressure at mesh interfaces (Pa) REAL, INTENT(IN) :: ppk(iip1, jjp1, llm) ! Exner at mid-layer REAL, INTENT(IN) :: flxw(iip1, jjp1, llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0) ! tendencies (in */s) from the physics REAL, INTENT(OUT) :: pdvfi(iip1, jjm, llm) ! tendency on covariant meridional wind REAL, INTENT(OUT) :: pdufi(iip1, jjp1, llm) ! tendency on covariant zonal wind REAL, INTENT(OUT) :: pdhfi(iip1, jjp1, llm) ! tendency on potential temperature (K/s) REAL, INTENT(OUT) :: pdqfi(iip1, jjp1, llm, nqtot) ! tendency on tracers REAL, INTENT(OUT) :: pdpsfi(iip1, jjp1) ! tendency on surface pressure (Pa/s) ! Local variables : ! ----------------- INTEGER :: i, j, l, ig0, ig, iq, itr REAL :: zpsrf(ngridmx) REAL :: zplev(ngridmx, llm + 1), zplay(ngridmx, llm) REAL :: zphi(ngridmx, llm), zphis(ngridmx) REAL :: zrot(iip1, jjm, llm) ! AdlC May 2014 REAL :: zufi(ngridmx, llm), zvfi(ngridmx, llm) REAL :: zrfi(ngridmx, llm) ! relative wind vorticity REAL :: ztfi(ngridmx, llm), zqfi(ngridmx, llm, nqtot) REAL :: zpk(ngridmx, llm) REAL :: pcvgu(ngridmx, llm), pcvgv(ngridmx, llm) REAL :: pcvgt(ngridmx, llm), pcvgq(ngridmx, llm, 2) REAL :: zdufi(ngridmx, llm), zdvfi(ngridmx, llm) REAL :: zdtfi(ngridmx, llm), zdqfi(ngridmx, llm, nqtot) REAL :: zdpsrf(ngridmx) REAL :: zdufic(ngridmx, llm), zdvfic(ngridmx, llm) REAL :: zdtfic(ngridmx, llm), zdqfic(ngridmx, llm, nqtot) REAL :: jH_cur_split, zdt_split LOGICAL :: debut_split, lafin_split INTEGER :: isplit REAL :: zsin(iim), zcos(iim), z1(iim) REAL :: zsinbis(iim), zcosbis(iim), z1bis(iim) REAL :: unskap, pksurcp REAL :: flxwfi(ngridmx, llm) ! Flux de masse verticale sur la grille physiq LOGICAL, SAVE :: firstcal = .TRUE., debut = .TRUE. ! REAL rdayvrai !----------------------------------------------------------------------- ! 1. Initialisations : ! -------------------- IF (firstcal) THEN debut = .TRUE. IF (ngridmx/=2 + (jjm - 1) * iim) THEN WRITE(lunout, *) 'STOP dans calfis' WRITE(lunout, *) & 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim' WRITE(lunout, *) ' ngridmx jjm iim ' WRITE(lunout, *) ngridmx, jjm, iim CALL abort_gcm("calfis", "", 1) ENDIF ELSE debut = .FALSE. ENDIF ! of IF (firstcal) !----------------------------------------------------------------------- ! 40. transformation des variables dynamiques en variables physiques: ! --------------------------------------------------------------- ! 41. pressions au sol (en Pascals) ! ---------------------------------- zpsrf(1) = pps(1, 1) ig0 = 2 DO j = 2, jjm CALL SCOPY(iim, pps(1, j), 1, zpsrf(ig0), 1) ig0 = ig0 + iim ENDDO zpsrf(ngridmx) = pps(1, jjp1) ! 42. pression intercouches et fonction d'Exner: ! ----------------------------------------------------------------- ! .... zplev definis aux (llm +1) interfaces des couches .... ! .... zplay definis aux ( llm ) milieux des couches .... ! ----------------------------------------------------------------- ! ... Exner = cp * ( p(l) / preff ) ** kappa .... unskap = 1. / kappa DO l = 1, llm zpk(1, l) = ppk(1, 1, l) zplev(1, l) = pp(1, 1, l) ig0 = 2 DO j = 2, jjm DO i = 1, iim zpk(ig0, l) = ppk(i, j, l) zplev(ig0, l) = pp(i, j, l) ig0 = ig0 + 1 ENDDO ENDDO zpk(ngridmx, l) = ppk(1, jjp1, l) zplev(ngridmx, l) = pp(1, jjp1, l) ENDDO zplev(1, llmp1) = pp(1, 1, llmp1) ig0 = 2 DO j = 2, jjm DO i = 1, iim zplev(ig0, llmp1) = pp(i, j, llmp1) ig0 = ig0 + 1 ENDDO ENDDO zplev(ngridmx, llmp1) = pp(1, jjp1, llmp1) ! ! 43. temperature naturelle (en K) et pressions milieux couches . ! --------------------------------------------------------------- DO l = 1, llm pksurcp = ppk(1, 1, l) / cpp zplay(1, l) = preff * pksurcp ** unskap ztfi(1, l) = pteta(1, 1, l) * pksurcp pcvgt(1, l) = pdteta(1, 1, l) * pksurcp / pmasse(1, 1, l) ig0 = 2 DO j = 2, jjm DO i = 1, iim pksurcp = ppk(i, j, l) / cpp zplay(ig0, l) = preff * pksurcp ** unskap ztfi(ig0, l) = pteta(i, j, l) * pksurcp pcvgt(ig0, l) = pdteta(i, j, l) * pksurcp / pmasse(i, j, l) ig0 = ig0 + 1 ENDDO ENDDO pksurcp = ppk(1, jjp1, l) / cpp zplay(ig0, l) = preff * pksurcp ** unskap ztfi (ig0, l) = pteta(1, jjp1, l) * pksurcp pcvgt(ig0, l) = pdteta(1, jjp1, l) * pksurcp / pmasse(1, jjp1, l) ENDDO ! 43.bis traceurs ! --------------- itr = 0 DO iq = 1, nqtot IF(.NOT.tracers(iq)%isAdvected) CYCLE itr = itr + 1 DO l = 1, llm zqfi(1, l, itr) = pq(1, 1, l, iq) ig0 = 2 DO j = 2, jjm DO i = 1, iim zqfi(ig0, l, itr) = pq(i, j, l, iq) ig0 = ig0 + 1 ENDDO ENDDO zqfi(ig0, l, itr) = pq(1, jjp1, l, iq) ENDDO ENDDO ! convergence dynamique pour les traceurs "EAU" ! Earth-specific treatment of first 2 tracers (water) IF (planet_type=="earth") THEN DO iq = 1, 2 DO l = 1, llm pcvgq(1, l, iq) = pdq(1, 1, l, iq) / pmasse(1, 1, l) ig0 = 2 DO j = 2, jjm DO i = 1, iim pcvgq(ig0, l, iq) = pdq(i, j, l, iq) / pmasse(i, j, l) ig0 = ig0 + 1 ENDDO ENDDO pcvgq(ig0, l, iq) = pdq(1, jjp1, l, iq) / pmasse(1, jjp1, l) ENDDO ENDDO endif ! of if (planet_type=="earth") ! Geopotentiel calcule par rapport a la surface locale: ! ----------------------------------------------------- CALL gr_dyn_fi(llm, iip1, jjp1, ngridmx, pphi, zphi) CALL gr_dyn_fi(1, iip1, jjp1, ngridmx, pphis, zphis) DO l = 1, llm DO ig = 1, ngridmx zphi(ig, l) = zphi(ig, l) - zphis(ig) ENDDO ENDDO ! .... Calcul de la vitesse verticale ( en Pa*m*s ou Kg/s ) .... ! JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux ! de masse est calclue dans advtrac.F ! DO l=1,llm ! pvervel(1,l)=pw(1,1,l) * g /apoln ! ig0=2 ! DO j=2,jjm ! DO i = 1, iim ! pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j) ! ig0 = ig0 + 1 ! ENDDO ! ENDDO ! pvervel(ig0,l)=pw(1,jjp1,l) * g /apols ! ENDDO ! 45. champ u: ! ------------ DO l = 1, llm DO j = 2, jjm ig0 = 1 + (j - 2) * iim zufi(ig0 + 1, l) = 0.5 * & (pucov(iim, j, l) / cu(iim, j) + pucov(1, j, l) / cu(1, j)) pcvgu(ig0 + 1, l) = 0.5 * & (pducov(iim, j, l) / cu(iim, j) + pducov(1, j, l) / cu(1, j)) DO i = 2, iim zufi(ig0 + i, l) = 0.5 * & (pucov(i - 1, j, l) / cu(i - 1, j) + pucov(i, j, l) / cu(i, j)) pcvgu(ig0 + i, l) = 0.5 * & (pducov(i - 1, j, l) / cu(i - 1, j) + pducov(i, j, l) / cu(i, j)) END DO END DO END DO ! Alvaro de la Camara (May 2014) ! 46.1 Calcul de la vorticite et passage sur la grille physique ! -------------------------------------------------------------- DO l = 1, llm DO i = 1, iim DO j = 1, jjm zrot(i, j, l) = (pvcov(i + 1, j, l) - pvcov(i, j, l) & + pucov(i, j + 1, l) - pucov(i, j, l)) & / (cu(i, j) + cu(i, j + 1)) & / (cv(i + 1, j) + cv(i, j)) * 4 enddo enddo ENDDO ! 46.champ v: ! ----------- DO l = 1, llm DO j = 2, jjm ig0 = 1 + (j - 2) * iim DO i = 1, iim zvfi(ig0 + i, l) = 0.5 * & (pvcov(i, j - 1, l) / cv(i, j - 1) + pvcov(i, j, l) / cv(i, j)) pcvgv(ig0 + i, l) = 0.5 * & (pdvcov(i, j - 1, l) / cv(i, j - 1) + pdvcov(i, j, l) / cv(i, j)) ENDDO zrfi(ig0 + 1, l) = 0.25 * (zrot(iim, j - 1, l) + zrot(iim, j, l) & + zrot(1, j - 1, l) + zrot(1, j, l)) DO i = 2, iim zrfi(ig0 + i, l) = 0.25 * (zrot(i - 1, j - 1, l) + zrot(i - 1, j, l) & + zrot(i, j - 1, l) + zrot(i, j, l)) ! AdlC MAY 2014 ENDDO ENDDO ENDDO ! 47. champs de vents aux pole nord ! ------------------------------ ! U = 1 / pi * integrale [ v * cos(long) * d long ] ! V = 1 / pi * integrale [ v * sin(long) * d long ] DO l = 1, llm z1(1) = (rlonu(1) - rlonu(iim) + 2. * pi) * pvcov(1, 1, l) / cv(1, 1) z1bis(1) = (rlonu(1) - rlonu(iim) + 2. * pi) * pdvcov(1, 1, l) / cv(1, 1) DO i = 2, iim z1(i) = (rlonu(i) - rlonu(i - 1)) * pvcov(i, 1, l) / cv(i, 1) z1bis(i) = (rlonu(i) - rlonu(i - 1)) * pdvcov(i, 1, l) / cv(i, 1) ENDDO DO i = 1, iim zcos(i) = COS(rlonv(i)) * z1(i) zcosbis(i) = COS(rlonv(i)) * z1bis(i) zsin(i) = SIN(rlonv(i)) * z1(i) zsinbis(i) = SIN(rlonv(i)) * z1bis(i) ENDDO zufi(1, l) = SSUM(iim, zcos, 1) / pi pcvgu(1, l) = SSUM(iim, zcosbis, 1) / pi zvfi(1, l) = SSUM(iim, zsin, 1) / pi pcvgv(1, l) = SSUM(iim, zsinbis, 1) / pi zrfi(1, l) = 0. ENDDO ! 48. champs de vents aux pole sud: ! --------------------------------- ! U = 1 / pi * integrale [ v * cos(long) * d long ] ! V = 1 / pi * integrale [ v * sin(long) * d long ] DO l = 1, llm z1(1) = (rlonu(1) - rlonu(iim) + 2. * pi) * pvcov(1, jjm, l) / cv(1, jjm) z1bis(1) = (rlonu(1) - rlonu(iim) + 2. * pi) * pdvcov(1, jjm, l) / cv(1, jjm) DO i = 2, iim z1(i) = (rlonu(i) - rlonu(i - 1)) * pvcov(i, jjm, l) / cv(i, jjm) z1bis(i) = (rlonu(i) - rlonu(i - 1)) * pdvcov(i, jjm, l) / cv(i, jjm) ENDDO DO i = 1, iim zcos(i) = COS(rlonv(i)) * z1(i) zcosbis(i) = COS(rlonv(i)) * z1bis(i) zsin(i) = SIN(rlonv(i)) * z1(i) zsinbis(i) = SIN(rlonv(i)) * z1bis(i) ENDDO zufi(ngridmx, l) = SSUM(iim, zcos, 1) / pi pcvgu(ngridmx, l) = SSUM(iim, zcosbis, 1) / pi zvfi(ngridmx, l) = SSUM(iim, zsin, 1) / pi pcvgv(ngridmx, l) = SSUM(iim, zsinbis, 1) / pi zrfi(ngridmx, l) = 0. ENDDO ! On change de grille, dynamique vers physiq, pour le flux de masse verticale CALL gr_dyn_fi(llm, iip1, jjp1, ngridmx, flxw, flxwfi) !----------------------------------------------------------------------- ! Appel de la physique: ! --------------------- ! WRITE(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys zdt_split = dtphys / nsplit_phys zdufic(:, :) = 0. zdvfic(:, :) = 0. zdtfic(:, :) = 0. zdqfic(:, :, :) = 0. IF (CPPKEY_PHYS) THEN DO isplit = 1, nsplit_phys jH_cur_split = jH_cur + (isplit - 1) * dtvr / (daysec * nsplit_phys) debut_split = debut.AND.isplit==1 lafin_split = lafin.AND.isplit==nsplit_phys ! if (planet_type=="earth") THEN CALL call_physiq(ngridmx, llm, nqtot, tracers(:)%name, & debut_split, lafin_split, & jD_cur, jH_cur_split, zdt_split, & zplev, zplay, & zpk, zphi, zphis, & presnivs, & zufi, zvfi, zrfi, ztfi, zqfi, & flxwfi, pducov, & zdufi, zdvfi, zdtfi, zdqfi, zdpsrf) ! ELSE IF ( planet_type=="generic" ) THEN ! CALL physiq (ngridmx, !! ngrid ! . llm, !! nlayer ! . nqtot, !! nq ! . tracers(:)%name,!! tracer names from dynamical core (given in infotrac) ! . debut_split, !! firstcall ! . lafin_split, !! lastcall ! . jD_cur, !! pday. see leapfrog ! . jH_cur_split, !! ptime "fraction of day" ! . zdt_split, !! ptimestep ! . zplev, !! pplev ! . zplay, !! pplay ! . zphi, !! pphi ! . zufi, !! pu ! . zvfi, !! pv ! . ztfi, !! pt ! . zqfi, !! pq ! . flxwfi, !! pw !! or 0. anyway this is for diagnostic. not used in physiq. ! . zdufi, !! pdu ! . zdvfi, !! pdv ! . zdtfi, !! pdt ! . zdqfi, !! pdq ! . zdpsrf, !! pdpsrf ! . tracerdyn) !! tracerdyn <-- utilite ??? ! ENDIF ! of if (planet_type=="earth") zufi(:, :) = zufi(:, :) + zdufi(:, :) * zdt_split zvfi(:, :) = zvfi(:, :) + zdvfi(:, :) * zdt_split ztfi(:, :) = ztfi(:, :) + zdtfi(:, :) * zdt_split zqfi(:, :, :) = zqfi(:, :, :) + zdqfi(:, :, :) * zdt_split zdufic(:, :) = zdufic(:, :) + zdufi(:, :) zdvfic(:, :) = zdvfic(:, :) + zdvfi(:, :) zdtfic(:, :) = zdtfic(:, :) + zdtfi(:, :) zdqfic(:, :, :) = zdqfic(:, :, :) + zdqfi(:, :, :) enddo ! of do isplit=1,nsplit_phys END IF zdufi(:, :) = zdufic(:, :) / nsplit_phys zdvfi(:, :) = zdvfic(:, :) / nsplit_phys zdtfi(:, :) = zdtfic(:, :) / nsplit_phys zdqfi(:, :, :) = zdqfic(:, :, :) / nsplit_phys !----------------------------------------------------------------------- ! transformation des tendances physiques en tendances dynamiques: ! --------------------------------------------------------------- ! tendance sur la pression : ! ----------------------------------- CALL gr_fi_dyn(1, ngridmx, iip1, jjp1, zdpsrf, pdpsfi) ! 62. enthalpie potentielle ! --------------------- DO l = 1, llm DO i = 1, iip1 pdhfi(i, 1, l) = cpp * zdtfi(1, l) / ppk(i, 1, l) pdhfi(i, jjp1, l) = cpp * zdtfi(ngridmx, l) / ppk(i, jjp1, l) ENDDO DO j = 2, jjm ig0 = 1 + (j - 2) * iim DO i = 1, iim pdhfi(i, j, l) = cpp * zdtfi(ig0 + i, l) / ppk(i, j, l) ENDDO pdhfi(iip1, j, l) = pdhfi(1, j, l) ENDDO ENDDO ! 62. humidite specifique ! --------------------- ! Ehouarn: removed this useless bit: was overwritten at step 63 anyways ! DO iq=1,nqtot ! DO l=1,llm ! DO i=1,iip1 ! pdqfi(i,1,l,iq) = zdqfi(1,l,iq) ! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq) ! ENDDO ! DO j=2,jjm ! ig0=1+(j-2)*iim ! DO i=1,iim ! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq) ! ENDDO ! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq) ! ENDDO ! ENDDO ! ENDDO ! 63. traceurs ! ------------ ! initialisation des tendances pdqfi(:, :, :, :) = 0. itr = 0 DO iq = 1, nqtot IF(.NOT.tracers(iq)%isAdvected) CYCLE itr = itr + 1 DO l = 1, llm DO i = 1, iip1 pdqfi(i, 1, l, iq) = zdqfi(1, l, itr) pdqfi(i, jjp1, l, iq) = zdqfi(ngridmx, l, itr) ENDDO DO j = 2, jjm ig0 = 1 + (j - 2) * iim DO i = 1, iim pdqfi(i, j, l, iq) = zdqfi(ig0 + i, l, itr) ENDDO pdqfi(iip1, j, l, iq) = pdqfi(1, j, l, itr) ENDDO ENDDO ENDDO ! 65. champ u: ! ------------ DO l = 1, llm DO i = 1, iip1 pdufi(i, 1, l) = 0. pdufi(i, jjp1, l) = 0. ENDDO DO j = 2, jjm ig0 = 1 + (j - 2) * iim DO i = 1, iim - 1 pdufi(i, j, l) = & 0.5 * (zdufi(ig0 + i, l) + zdufi(ig0 + i + 1, l)) * cu(i, j) ENDDO pdufi(iim, j, l) = & 0.5 * (zdufi(ig0 + 1, l) + zdufi(ig0 + iim, l)) * cu(iim, j) pdufi(iip1, j, l) = pdufi(1, j, l) ENDDO ENDDO ! 67. champ v: ! ------------ DO l = 1, llm DO j = 2, jjm - 1 ig0 = 1 + (j - 2) * iim DO i = 1, iim pdvfi(i, j, l) = & 0.5 * (zdvfi(ig0 + i, l) + zdvfi(ig0 + i + iim, l)) * cv(i, j) ENDDO pdvfi(iip1, j, l) = pdvfi(1, j, l) ENDDO ENDDO ! 68. champ v pres des poles: ! --------------------------- ! v = U * cos(long) + V * SIN(long) DO l = 1, llm DO i = 1, iim pdvfi(i, 1, l) = & zdufi(1, l) * COS(rlonv(i)) + zdvfi(1, l) * SIN(rlonv(i)) pdvfi(i, jjm, l) = zdufi(ngridmx, l) * COS(rlonv(i)) & + zdvfi(ngridmx, l) * SIN(rlonv(i)) pdvfi(i, 1, l) = & 0.5 * (pdvfi(i, 1, l) + zdvfi(i + 1, l)) * cv(i, 1) pdvfi(i, jjm, l) = & 0.5 * (pdvfi(i, jjm, l) + zdvfi(ngridmx - iip1 + i, l)) * cv(i, jjm) ENDDO pdvfi(iip1, 1, l) = pdvfi(1, 1, l) pdvfi(iip1, jjm, l) = pdvfi(1, jjm, l) ENDDO !----------------------------------------------------------------------- firstcal = .FALSE. RETURN END SUBROUTINE calfis