Index: LMDZ.3.3/branches/rel-LF/gcm.def =================================================================== --- LMDZ.3.3/branches/rel-LF/gcm.def (revision 265) +++ LMDZ.3.3/branches/rel-LF/gcm.def (revision 265) @@ -0,0 +1,81 @@ +# +## $Header$ +# +## Jour de l'etat initial ( = 350 si 20 Decembre ,par expl. ,comme ici ) +dayref=1 +## Annee de l'etat initial ( avec 4 chiffres ) +anneeref=1979 +## Nombre de jours d'integration +nday=5 +## nombre de pas par jour (multiple de iperiod) ( ici pour dt = 1 min ) +day_step=480 +## periode pour le pas Matsuno (en pas) +iperiod=5 +## periode de sortie des variables de controle (en pas) +iconser=240 +## periode d'ecriture du fichier histoire (en jour) +iecri=1 +## periode de stockage fichier histmoy (en jour) +periodav=1. +## periode de la dissipation (en pas) +idissip=5 +## choix de l'operateur de dissipation (star ou non star ) +lstardis=y +## nombre d'iterations de l'operateur de dissipation gradiv +nitergdiv=1 +## nombre d'iterations de l'operateur de dissipation nxgradrot +nitergrot=2 +## nombre d'iterations de l'operateur de dissipation divgrad +niterh=2 +## temps de dissipation des plus petites long.d ondes pour u,v (gradiv) +tetagdiv=3600. +## temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot) +tetagrot=3600. +## temps de dissipation des plus petites long.d ondes pour h ( divgrad) +tetatemp=3600. +## coefficient pour gamdissip +coefdis=0. +## choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) +purmats=n +## avec ou sans physique +physic=y +## periode de la physique (en pas) +iphysiq=10 +## frequence (en jours ) de l'ecriture du fichier histphy +ecritphy=30 +## Cycle diurne ou non +cycle_diurne=y +## Soil Model ou non +soil_model=y +## Choix ou non de New oliq +new_oliq=y +## Orodr ou non pour l orographie +ok_orodr=y +## Orolf ou non pour l orographie +ok_orolf=y +## Si = .T. , lecture du fichier limit avec la bonne annee +ok_limitvrai=n +## Nombre d'appels des routines de rayonnements ( par jour) +nbapp_rad=12 +## Flag pour la convection ( 1 pour LMD, 2 pour Tiedtke, 3 CCM ) +iflag_con=2 +## longitude en degres du centre du zoom +clon=0. +## latitude en degres du centre du zoom +clat=0. +## facteur de grossissement du zoom,selon longitude +grossismx=1.0 +## facteur de grossissement du zoom ,selon latitude +grossismy=1.0 +## Fonction f(y) hyperbolique si = .true. , sinon sinusoidale +fxyhypb=y +## extension en longitude de la zone du zoom ( fraction de la zone totale) +dzoomx=0.0 +## extension en latitude de la zone du zoom ( fraction de la zone totale) +dzoomy=0.0 +##raideur du zoom en X +taux=3. +##raideur du zoom en Y +tauy=3. +## Fonction f(y) avec y = Sin(latit.) si = .true. , sinon y = latit. +ysinus=y Index: LMDZ.3.3/branches/rel-LF/libf/bibio/writephys.F90 =================================================================== --- LMDZ.3.3/branches/rel-LF/libf/bibio/writephys.F90 (revision 265) +++ LMDZ.3.3/branches/rel-LF/libf/bibio/writephys.F90 (revision 265) @@ -0,0 +1,369 @@ +! +! $Header$ +! + MODULE writephys + +! +! Wrapper d'IOIPSL pour l'ecriture des sorties de la physique. Un seul appel +! par variable devrait (?) suffire pour l'initialisation et l'ecriture +! +! LF, LMD, 07/2001 +! +! 2 routines principales: +! writephy_ini pour l'initialisation des fichiers (appel histbeg, histvert) +! writephy pour l'ecriture des champs (appel histdef, histwrite) +! dans le futur: +! writephy_def pour definir des profils de variables +! + USE IOIPSL + + IMPLICIT none + + PRIVATE + PUBLIC :: writephy_ini, writephy, writephy_def, writephy_sync + +! +! variables locales +! +! nombre de fichiers maximum a ouvrir + integer, parameter :: nb_files_max = 20 +! nombre de variables maximum par fichier + integer, parameter :: nb_var_max = 500 +! nombre maximum de "profils" de variables + integer, parameter :: nb_prof_max = 10 +! structure d'infos sur les fichiers + type fichier + character*30 :: file_name + integer :: nitau + integer :: file_id + integer :: isize, jsize, lsize, phy_lon + integer :: nhori, nvert + logical :: define + end type fichier + type (fichier), dimension(nb_files_max), save :: hist_files +! structure de profils de variables + type profils + real :: freq_op + real :: freq_wri + character*6 :: var_op + integer :: zsize + integer :: reg_size + integer,dimension(:), pointer :: reg_index + end type profils + type (profils), dimension(nb_prof_max), save :: var_profs + +! liste des variables par fichier + character (len=10), save, dimension(nb_var_max,nb_files_max) :: list_var='undefined' +! nombre de variables par fichier + integer, save, dimension(nb_files_max) :: nb_var = 0 +! nombre de bits à sauvegarder dans les fichiers + integer, save :: nbits = 32 +! +! Quelques initialisations +! + + CONTAINS + +! +!########################################################################### +! + SUBROUTINE writephy_ini(nid_file, nom_fichier, klon, iim, jjm, llm, & + & rlon, rlat, zlev, & + & date0, dtime) +! +! Initialisation des fichiers histoire de la physique +! Appels a histbeg, histvert +! Remplissage des structures d'informations sur les fichiers +! +! Entree: +! nid_file numero/index du fichier a initialiser +! nom_fichier nom du fichier +! klon taille des champs physiques +! iim, jjm, llm taille en i,j,k +! rlon, rlat, zlev coordonnees en x, y, z +! date0 date debut +! dtime pas de temps elementaire de la physique +! +! Declaration des parametres d'entree + integer :: nid_file + character (len=30) :: nom_fichier + integer :: iim, jjm, llm, klon + real, dimension(iim) :: rlon + real, dimension(jjm) :: rlat + real, dimension(llm) :: zlev + real :: date0, dtime +! +! Variables locales +! + integer :: i, file_id, nvert, nhori + real, dimension(iim, jjm) :: temp_lon, temp_lat + logical, save :: first_call = .true. + integer,parameter :: nulou = 6 +! +! Quelques initialisations +! + if (first_call) then + hist_files(:)%define = .true. + var_profs(:)%freq_wri = -999 + hist_files(:)%nitau = 1 + first_call= .false. + endif +! +! Mise des coordonnees sur une grille 2D +! + call gr_fi_ecrit(1,klon,iim,jjm,rlon,temp_lon) + do i = 1, iim + temp_lon(i,1) = rlon(i+1) + temp_lon(i,jjm) = rlon(i+1) + enddo + call gr_fi_ecrit(1,klon,iim,jjm,rlat,temp_lat) +! +! Initialisation du fichier +! + call histbeg(nom_fichier, iim, temp_lon, jjm, temp_lat, & + & 1, iim, 1, jjm, 0, date0, dtime, & + & nhori, file_id) + call histvert(file_id, "presnivs", "Vertical levels", "mb", & + & llm, zlev, nvert) + +! +! Remplissage des infos +! + hist_files(nid_file)%file_name = nom_fichier + hist_files(nid_file)%file_id = file_id + hist_files(nid_file)%isize = iim + hist_files(nid_file)%jsize = jjm + hist_files(nid_file)%lsize = llm + hist_files(nid_file)%phy_lon = klon + hist_files(nid_file)%nhori = nhori + hist_files(nid_file)%nvert = nvert + + write(nulou,*)'####################################' + write(nulou,*)' Definition du fichier no ',nid_file + write(nulou,*)' nom du fichier ',hist_files(nid_file)%file_name + write(nulou,*)' identifiant ',hist_files(nid_file)%file_id + write(nulou,*)' taille en i ',hist_files(nid_file)%isize + write(nulou,*)' taille en j ',hist_files(nid_file)%jsize + write(nulou,*)' taille en l ',hist_files(nid_file)%lsize + +! + END SUBROUTINE writephy_ini +! +!########################################################################### +! + SUBROUTINE writephy_def(profil_n, var_op, freq_op, freq_wri, zsize, & + & reg_size, reg_index) +! +! Definition de profil type de sortie de variables +! (moyenne, instantanee, indexee ...) +! +! Parametres d'entree +! profil_n numero du profil-type +! var_op operation a effectuer sur la variable +! freq_op frequence de l'operation +! freq_wri frequence d'ecriture +! zsize variable 2D/3D +! reg_size taille de la region a sortir +! reg_index indices de la region a sortir +! + integer :: profil_n + character (len = 6) :: var_op + real :: freq_op, freq_wri + integer :: zsize + integer, optional :: reg_size + integer, dimension(*), optional :: reg_index +! +! variables locales +! + character (len = 20) :: modname = 'writephy_def' + character (len = 80) :: message + integer :: dummy_size + integer :: error + integer,parameter :: nulou = 6 +! +! Differentes verif +! + if (profil_n > nb_prof_max) then + message = 'numero de profil de variable > nbre maximal de profil' + call abort_gcm(modname, message, 1) + endif + if (var_profs(profil_n)%freq_wri /= -999) then + message = 'numero de profil deja attribue' + call abort_gcm(modname, message, 1) + endif +! +! Remplissage structure infos +! + var_profs(profil_n)%var_op = var_op + var_profs(profil_n)%freq_op = freq_op + var_profs(profil_n)%freq_wri = freq_wri + var_profs(profil_n)%zsize = zsize +! +! test pour region +! + if (present (reg_size)) then + if ( .not. present (reg_index)) then + message = 'reg_size defini mais pas de region definie' + call abort_gcm(modname, message, 1) + endif + allocate(var_profs(profil_n)%reg_index(reg_size), stat = error) + if (error /= 0) then + message='Pb allocation reg_index' + call abort_gcm(modname,message,1) + endif + var_profs(profil_n)%reg_size = reg_size + var_profs(profil_n)%reg_index = reg_index(1:reg_size) + else + dummy_size = 1 + allocate(var_profs(profil_n)%reg_index(dummy_size), stat = error) + if (error /= 0) then + message='Pb allocation reg_index' + call abort_gcm(modname,message,1) + endif + var_profs(profil_n)%reg_size = dummy_size + var_profs(profil_n)%reg_index = 0 + endif + + write(nulou,*)' Definition du profil de variable numero ', profil_n + write(nulou,*)' operation ',var_profs(profil_n)%var_op + write(nulou,*)' frequence d''operation ',var_profs(profil_n)%freq_op + write(nulou,*)' frequence d''ecriture ',var_profs(profil_n)%freq_wri + write(nulou,*)' 2D/3D ',var_profs(profil_n)%zsize + write(nulou,*)' taille de la region ',var_profs(profil_n)%reg_size + + END SUBROUTINE writephy_def +! +!########################################################################### +! +! A faire: rendre var_title et var_units optionels par lecture d'un tableau +! +! + SUBROUTINE writephy(file, iprof, var_name, data, var_title, var_units) +! +! Definition et ecriture des variables +! +! file numero du fichier dans lequel ecrire +! iprof profil de la variable a ecrire +! var_name nom de la variable a ecrire +! data les donnees effectives a ecrire +! var_title "vrai" nom de la variable, optionel, si non present +! interrogation d'un tableau +! var_units unite de la variable, optionel, si non present +! interrogation d'un tableau + + integer :: file, iprof + character (len=10) :: var_name + real, dimension(*) :: data + character (len=40) :: var_title + character (len=20) :: var_units +! +! variables locales +! + integer :: i, error + character (len=6) :: var_op + real :: freq_op, freq_wri + integer :: file_id, isize, jsize, lsize, phy_lon, zsize + integer :: nhori, nvert, klon + integer :: nitau + real, dimension(:,:,:), allocatable :: temp_data + character (len = 20) :: modname = 'writephy' + character (len = 80) :: message +! +! Test: toujours en mode definition des variables? +! + if (var_name == list_var(1, file)) then + hist_files(file)%nitau = hist_files(file)%nitau + 1 + endif + if (hist_files(file)%define) then + if (var_name == list_var(1, file)) then +! on a fait le tour des variables + hist_files(file)%define = .false. + call histend(hist_files(file)%file_id) +! hist_files(file)%nitau = hist_files(file)%nitau + 1 + else +! +! pour que l'appel a histdef soit plus lisible, on range tout dans des tampons +! + file_id = hist_files(file)%file_id + isize = hist_files(file)%isize + jsize = hist_files(file)%jsize + nhori = hist_files(file)%nhori + klon = hist_files(file)%phy_lon + var_op = var_profs(iprof)%var_op + freq_op = var_profs(iprof)%freq_op + freq_wri = var_profs(iprof)%freq_wri + if (var_profs(iprof)%zsize == 0) then + lsize = 1 + nvert = -99 + else + lsize = hist_files(file)%lsize + nvert = hist_files(file)%nvert + endif + call histdef(file_id, var_name, var_title, var_units, & + & isize, jsize, nhori, lsize, 1, lsize, nvert, nbits, & + & var_op, freq_op, freq_wri) + nb_var(file) = nb_var(file) + 1 + list_var(nb_var(file), file) = var_name + endif + endif +! +! On ecrit la variable +! +! Preparation +! + file_id = hist_files(file)%file_id + isize = hist_files(file)%isize + jsize = hist_files(file)%jsize + nhori = hist_files(file)%nhori + var_op = var_profs(iprof)%var_op + freq_op = var_profs(iprof)%freq_op + freq_wri = var_profs(iprof)%freq_wri + nitau = hist_files(file)%nitau + + if (var_profs(iprof)%zsize == 0) then + lsize = 1 + else + lsize = hist_files(file)%lsize + endif + allocate(temp_data(isize,jsize,lsize), stat = error) + if (error /= 0) then + message='Pb allocation temp_data' + call abort_gcm(modname, message, 1) + endif +! +! Ecriture +! + call gr_fi_ecrit(lsize, klon, isize, jsize, data, temp_data) + call histwrite(file_id, var_name,nitau,temp_data, & + & var_profs(iprof)%reg_size,var_profs(iprof)%reg_index) + + return +! + END SUBROUTINE writephy +! +!########################################################################### +! + SUBROUTINE writephy_sync(nid_file) +! +! Flush des donnees dans le fichier et eventuellement fermeture du +! mode 'define' +! +! Entree: +! nid_file numero du fichier a traiter +! + integer :: nid_file +! + if (hist_files(nid_file)%define) then + call histend(hist_files(nid_file)%file_id) + hist_files(nid_file)%define = .false. + endif + + call histsync(hist_files(nid_file)%file_id) + +return + + END SUBROUTINE writephy_sync +! +!########################################################################### +! + END MODULE writephys Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/conf_dat3d.F =================================================================== --- LMDZ.3.3/branches/rel-LF/libf/dyn3d/conf_dat3d.F (revision 265) +++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/conf_dat3d.F (revision 265) @@ -0,0 +1,293 @@ +C +C $Header$ +C + SUBROUTINE conf_dat3d( title, lons,lats,levs,xd,yd,zd,xf,yf,zf, + , champd , interbar ) +c +c Auteur : P. Le Van +c +c Ce s-pr. configure le champ de donnees 3D 'champd' de telle facon +c qu'on ait - pi a pi en longitude +c qu'on ait pi/2. a - pi/2. en latitude +c et qu'on ait les niveaux verticaux variant du sol vers le ht de l'atmos. +c ( en Pascals ) . +c +c xd et yd sont les longitudes et latitudes initiales +c zd les pressions initiales +c +c xf et yf sont les longitudes et latitudes en sortie , eventuellement +c modifiees pour etre configurees comme ci-dessus . +c zf les pressions en sortie +c +c champd en meme temps le champ initial et final +c +c interbar = .TRUE. si on appelle l'interpo. barycentrique inter_barxy +c sinon , l'interpolation grille_m ( grid_atob ) . +c + + IMPLICIT NONE + +c *** Arguments en entree *** + CHARACTER*(*) :: title + INTEGER lons, lats, levs + REAL xd(lons), yd(lats), zd(levs) + LOGICAL interbar +c +c *** Arguments en sortie *** + REAL xf(lons), yf(lats), zf(levs) + +c *** Arguments en entree et sortie *** + REAL champd(lons,lats,levs) + +c *** Variables locales *** +c + REAL pi,pis2,depi,presmax + LOGICAL radianlon, invlon ,radianlat, invlat, invlev, alloc + REAL rlatmin,rlatmax,oldxd1 + INTEGER i,j,ip180,ind,l + + REAL, ALLOCATABLE :: xtemp(:) + REAL, ALLOCATABLE :: ytemp(:) + REAL, ALLOCATABLE :: ztemp(:) + REAL, ALLOCATABLE :: champf(:,:,:) + + +c WRITE(6,*) ' Conf_dat3d pour ',title + + ALLOCATE(xtemp(lons)) + ALLOCATE(ytemp(lats)) + ALLOCATE(ztemp(levs)) + + DO i = 1, lons + xtemp(i) = xd(i) + ENDDO + DO j = 1, lats + ytemp(j) = yd(j) + ENDDO + DO l = 1, levs + ztemp(l) = zd(l) + ENDDO + + pi = 2. * ASIN(1.) + pis2 = pi/2. + depi = 2. * pi + + IF( xtemp(1).GE.-pi-0.5.AND. xtemp(lons).LE.pi+0.5 ) THEN + radianlon = .TRUE. + invlon = .FALSE. + ELSE IF (xtemp(1).GE.-0.5.AND.xtemp(lons).LE.depi+0.5 ) THEN + radianlon = .TRUE. + invlon = .TRUE. + ELSE IF ( xtemp(1).GE.-180.5.AND. xtemp(lons).LE.180.5 ) THEN + radianlon = .FALSE. + invlon = .FALSE. + ELSE IF ( xtemp(1).GE.-0.5.AND.xtemp(lons).LE.360.5 ) THEN + radianlon = .FALSE. + invlon = .TRUE. + ELSE + WRITE(6,*) 'Pbs. sur les longitudes des donnees pour le fichier' + , , title + ENDIF + + invlat = .FALSE. + + IF( ytemp(1).LT.ytemp(lats) ) THEN + invlat = .TRUE. + ENDIF + + rlatmin = MIN( ytemp(1), ytemp(lats) ) + rlatmax = MAX( ytemp(1), ytemp(lats) ) + + IF( rlatmin.GE.-pis2-0.5.AND.rlatmax.LE.pis2+0.5)THEN + radianlat = .TRUE. + ELSE IF ( rlatmin.GE.-90.-0.5.AND.rlatmax.LE.90.+0.5 ) THEN + radianlat = .FALSE. + ELSE + WRITE(6,*) ' Pbs. sur les latitudes des donnees pour le fichier' + , , title + ENDIF + + IF( .NOT. radianlon ) THEN + DO i = 1, lons + xtemp(i) = xtemp(i) * pi/180. + ENDDO + ENDIF + + IF( .NOT. radianlat ) THEN + DO j = 1, lats + ytemp(j) = ytemp(j) * pi/180. + ENDDO + ENDIF + + + alloc =.FALSE. + + IF ( invlon ) THEN + + ALLOCATE(champf(lons,lats,levs)) + alloc = .TRUE. + + DO i = 1 ,lons + xf(i) = xtemp(i) + ENDDO + + DO l = 1, levs + DO j = 1, lats + DO i= 1, lons + champf (i,j,l) = champd (i,j,l) + ENDDO + ENDDO + ENDDO +c +c *** On tourne les longit. pour avoir - pi a + pi **** +c + DO i=1,lons + IF( xf(i).GT. pi ) THEN + GO TO 88 + ENDIF + ENDDO + +88 CONTINUE +c + ip180 = i + + DO i = 1,lons + IF (xf(i).GT. pi) THEN + xf(i) = xf(i) - depi + ENDIF + ENDDO + + DO i= ip180,lons + ind = i-ip180 +1 + xtemp(ind) = xf(i) + ENDDO + + DO i= ind +1,lons + xtemp(i) = xf(i-ind) + ENDDO + +c ..... on tourne les longitudes pour champf .... +c + DO l = 1,levs + DO j = 1,lats + DO i = ip180,lons + ind = i-ip180 +1 + champd (ind,j,l) = champf (i,j,l) + ENDDO + + DO i= ind +1,lons + champd (i,j,l) = champf (i-ind,j,l) + ENDDO + ENDDO + ENDDO + + DEALLOCATE(xtemp) + ENDIF +c +c ***** fin de IF(invlon) **** + + IF ( invlat ) THEN + + IF(.NOT.alloc) THEN + ALLOCATE(champf(lons,lats,levs)) + alloc = .TRUE. + ENDIF + + DO j = 1, lats + yf(j) = ytemp(j) + ENDDO + + DO l = 1,levs + DO j = 1, lats + DO i = 1,lons + champf(i,j,l) = champd(i,j,l) + ENDDO + ENDDO + + DO j = 1, lats + ytemp( lats-j+1 ) = yf(j) + DO i = 1, lons + champd (i,lats-j+1,l) = champf (i,j,l) + ENDDO + ENDDO + ENDDO + + DEALLOCATE(ytemp) + + ENDIF + +c ***** fin de IF(invlat) **** +c +c + IF( interbar ) THEN + oldxd1 = xtemp(1) + DO i = 1, lons -1 + xtemp(i) = 0.5 * ( xtemp(i) + xtemp(i+1) ) + ENDDO + xtemp(lons) = 0.5 * ( xtemp(lons) + oldxd1 + depi ) + + DO j = 1, lats -1 + ytemp(j) = 0.5 * ( ytemp(j) + ytemp(j+1) ) + ENDDO + ENDIF +c + + invlev = .FALSE. + IF( ztemp(1).LT.ztemp(levs) ) invlev = .TRUE. + + presmax = MAX( ztemp(1), ztemp(levs) ) + IF( presmax.LT.1200. ) THEN + DO l = 1,levs + ztemp(l) = ztemp(l) * 100. + ENDDO + ENDIF + + IF( invlev ) THEN + + IF(.NOT.alloc) THEN + ALLOCATE(champf(lons,lats,levs)) + alloc = .TRUE. + ENDIF + + DO l = 1,levs + zf(l) = ztemp(l) + ENDDO + + DO l = 1,levs + DO j = 1, lats + DO i = 1,lons + champf(i,j,l) = champd(i,j,l) + ENDDO + ENDDO + ENDDO + + DO l = 1,levs + ztemp(levs+1-l) = zf(l) + ENDDO + + DO l = 1,levs + DO j = 1, lats + DO i = 1,lons + champd(i,j,levs+1-l) = champf(i,j,l) + ENDDO + ENDDO + ENDDO + + DEALLOCATE(ztemp) + + ENDIF + + IF(alloc) DEALLOCATE(champf) + + DO i = 1, lons + xf(i) = xtemp(i) + ENDDO + DO j = 1, lats + yf(j) = ytemp(j) + ENDDO + DO l = 1, levs + zf(l) = ztemp(l) + ENDDO + + RETURN + END Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/extrapol.F =================================================================== --- LMDZ.3.3/branches/rel-LF/libf/dyn3d/extrapol.F (revision 265) +++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/extrapol.F (revision 265) @@ -0,0 +1,198 @@ +C +C $Header$ +C + SUBROUTINE extrapol (pfild, kxlon, kylat, pmask, + . norsud, ldper, knbor, pwork) + IMPLICIT none +c +c OASIS routine (Adaptation: Laurent Li, le 14 mars 1997) +c Fill up missed values by using the neighbor points +c + INTEGER kxlon, kylat ! longitude and latitude dimensions (Input) + INTEGER knbor ! minimum neighbor number (Input) + LOGICAL norsud ! True if field is from North to South (Input) + LOGICAL ldper ! True if take into account the periodicity (Input) + REAL pmask ! mask value (Input) + REAL pfild(kxlon,kylat) ! field to be extrapolated (Input/Output) + REAL pwork(kxlon,kylat) ! working space +c + REAL zwmsk + INTEGER incre, idoit, i, j, k, inbor, ideb, ifin, ilon, jlat + INTEGER ix(9), jy(9) ! index arrays for the neighbors coordinates + REAL zmask(9) +C +C We search over the eight closest neighbors +C +C j+1 7 8 9 +C j 4 5 6 Current point 5 --> (i,j) +C j-1 1 2 3 +C i-1 i i+1 +c +c + IF (norsud) THEN + DO j = 1, kylat + DO i = 1, kxlon + pwork(i,j) = pfild(i,kylat-j+1) + ENDDO + ENDDO + DO j = 1, kylat + DO i = 1, kxlon + pfild(i,j) = pwork(i,j) + ENDDO + ENDDO + ENDIF +c + incre = 0 +c + DO j = 1, kylat + DO i = 1, kxlon + pwork(i,j) = pfild(i,j) + ENDDO + ENDDO +c +C* To avoid problems in floating point tests + zwmsk = pmask - 1.0 +c +200 CONTINUE + incre = incre + 1 + DO 99999 j = 1, kylat + DO 99999 i = 1, kxlon + IF (pfild(i,j).GT. zwmsk) THEN + pwork(i,j) = pfild(i,j) + inbor = 0 + ideb = 1 + ifin = 9 +C +C* Fill up ix array + ix(1) = MAX (1,i-1) + ix(2) = i + ix(3) = MIN (kxlon,i+1) + ix(4) = MAX (1,i-1) + ix(5) = i + ix(6) = MIN (kxlon,i+1) + ix(7) = MAX (1,i-1) + ix(8) = i + ix(9) = MIN (kxlon,i+1) +C +C* Fill up iy array + jy(1) = MAX (1,j-1) + jy(2) = MAX (1,j-1) + jy(3) = MAX (1,j-1) + jy(4) = j + jy(5) = j + jy(6) = j + jy(7) = MIN (kylat,j+1) + jy(8) = MIN (kylat,j+1) + jy(9) = MIN (kylat,j+1) +C +C* Correct latitude bounds if southernmost or northernmost points + IF (j .EQ. 1) ideb = 4 + IF (j .EQ. kylat) ifin = 6 +C +C* Account for periodicity in longitude +C + IF (ldper) THEN + IF (i .EQ. kxlon) THEN + ix(3) = 1 + ix(6) = 1 + ix(9) = 1 + ELSE IF (i .EQ. 1) THEN + ix(1) = kxlon + ix(4) = kxlon + ix(7) = kxlon + ENDIF + ELSE + IF (i .EQ. 1) THEN + ix(1) = i + ix(2) = i + 1 + ix(3) = i + ix(4) = i + 1 + ix(5) = i + ix(6) = i + 1 + ENDIF + IF (i .EQ. kxlon) THEN + ix(1) = i -1 + ix(2) = i + ix(3) = i - 1 + ix(4) = i + ix(5) = i - 1 + ix(6) = i + ENDIF +C + IF (i .EQ. 1 .OR. i .EQ. kxlon) THEN + jy(1) = MAX (1,j-1) + jy(2) = MAX (1,j-1) + jy(3) = j + jy(4) = j + jy(5) = MIN (kylat,j+1) + jy(6) = MIN (kylat,j+1) +C + ideb = 1 + ifin = 6 + IF (j .EQ. 1) ideb = 3 + IF (j .EQ. kylat) ifin = 4 + ENDIF + ENDIF ! end for ldper test +C +C* Find unmasked neighbors +C + DO 230 k = ideb, ifin + zmask(k) = 0. + ilon = ix(k) + jlat = jy(k) + IF (pfild(ilon,jlat) .LT. zwmsk) THEN + zmask(k) = 1. + inbor = inbor + 1 + ENDIF + 230 CONTINUE +C +C* Not enough points around point P are unmasked; interpolation on P +C will be done in a future call to extrap. +C + IF (inbor .GE. knbor) THEN + pwork(i,j) = 0. + DO k = ideb, ifin + ilon = ix(k) + jlat = jy(k) + pwork(i,j) = pwork(i,j) + $ + pfild(ilon,jlat) * zmask(k)/FLOAT(inbor) + ENDDO + ENDIF +C + ENDIF +99999 CONTINUE +C +C* 3. Writing back unmasked field in pfild +C ------------------------------------ +C +C* pfild then contains: +C - Values which were not masked +C - Interpolated values from the inbor neighbors +C - Values which are not yet interpolated +C + idoit = 0 + DO j = 1, kylat + DO i = 1, kxlon + IF (pwork(i,j) .GT. zwmsk) idoit = idoit + 1 + pfild(i,j) = pwork(i,j) + ENDDO + ENDDO +c + IF (idoit .ne. 0) GOTO 200 +ccc PRINT*, "Number of extrapolation steps incre =", incre +c + IF (norsud) THEN + DO j = 1, kylat + DO i = 1, kxlon + pwork(i,j) = pfild(i,kylat-j+1) + ENDDO + ENDDO + DO j = 1, kylat + DO i = 1, kxlon + pfild(i,j) = pwork(i,j) + ENDDO + ENDDO + ENDIF +c + RETURN + END Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/limit_netcdf.F =================================================================== --- LMDZ.3.3/branches/rel-LF/libf/dyn3d/limit_netcdf.F (revision 265) +++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/limit_netcdf.F (revision 265) @@ -0,0 +1,1062 @@ +C +C $Header$ +C + SUBROUTINE limit_netcdf ( interbar, extrap, oldice ) +c + IMPLICIT none +c +c------------------------------------------------------------- +C Author : L. Fairhead +C Date : 27/01/94 +C Objet : Construction des fichiers de conditions aux limites +C pour le nouveau +C modele a partir de fichiers de climatologie. Les deux +C grilles doivent etre regulieres +c +c Modifie par z.x.li (le23mars1994) +c Modifie par L. Fairhead (fairhead@lmd.jussieu.fr) septembre 1999 +c pour lecture netcdf dans LMDZ.3.3 +c Modifie par P;Le Van , juillet 2001 +c------------------------------------------------------------- +c +#include "dimensions.h" +#include "paramet.h" +#include "control.h" +#include "logic.h" +#include "netcdf.inc" +#include "comvert.h" +#include "comgeom2.h" +#include "comconst.h" +c +c----------------------------------------------------------------------- + LOGICAL interbar, extrap, oldice + + INTEGER KIDIA, KFDIA, KLON, KLEV + PARAMETER (KIDIA=1,KFDIA=iim*(jjm-1)+2, + . KLON=KFDIA-KIDIA+1,KLEV=llm) +c----------------------------------------------------------------------- + REAL phy_nat(klon,360), phy_nat0(klon) + REAL phy_alb(klon,360) + REAL phy_sst(klon,360) + REAL phy_bil(klon,360) + REAL phy_rug(klon,360) + REAL phy_ice(klon,360) +c + REAL masque(iip1,jjp1) + REAL mask(iim,jjp1) + +C Declarations pour le champ de depart + INTEGER imdep, jmdep,lmdep + INTEGER tbid + PARAMETER ( tbid = 60 ) ! >52 semaines + REAL timecoord(tbid) +c + REAL , ALLOCATABLE :: dlon_msk(:), dlat_msk(:) + REAL , ALLOCATABLE :: lonmsk_ini(:), latmsk_ini(:) + REAL , ALLOCATABLE :: dlon(:), dlat(:) + REAL , ALLOCATABLE :: dlon_ini(:), dlat_ini(:) + REAL , ALLOCATABLE :: champ_msk(:), champ(:) + REAL , ALLOCATABLE :: work(:,:) + + CHARACTER*25 title + +C Declarations pour le champ interpole 2D + REAL champint(iim,jjp1) + real chmin,chmax + +C Declarations pour le champ interpole 3D + REAL champtime(iim,jjp1,tbid) + REAL timeyear(tbid) + REAL champan(iip1,jjp1,366) + +C Declarations pour l'inteprolation verticale + REAL ax(tbid), ay(tbid) + REAL by + REAL yder(tbid) + + + INTEGER ierr + INTEGER dimfirst(3) + INTEGER dimlast(3) +c + INTEGER nid, ndim, ntim + INTEGER dims(2), debut(2), epais(2) + INTEGER id_tim + INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB + + INTEGER i, j, k, l +c Diverses variables locales + REAL time + + INTEGER longcles + PARAMETER ( longcles = 20 ) + REAL clesphy0 ( longcles ) +#include "serre.h" + INTEGER ncid,varid,ndimid(4),dimid + character*30 namedim + + pi = 4. * ATAN(1.) + rad = 6 371 229. + omeg = 4.* ASIN(1.)/(24.*3600.) + g = 9.8 + daysec = 86400. + kappa = 0.2857143 + cpp = 1004.70885 + dtvr = daysec/FLOAT(day_step) +c +C Traitement du relief au sol +c + write(*,*) 'Traitement du relief au sol pour fabriquer masque' + ierr = NF_OPEN('Relief.nc', NF_NOWRITE, ncid) + + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ierr = NF_INQ_VARID(ncid,'RELIEF',varid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_VARDIMID (ncid,varid,ndimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, 'dimension ', imdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ALLOCATE( lonmsk_ini(imdep) ) + ALLOCATE( dlon_msk(imdep) ) + +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,lonmsk_ini) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,lonmsk_ini) +#endif + +c + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, 'dimension ', jmdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ALLOCATE( latmsk_ini(jmdep) ) + ALLOCATE( dlat_msk(jmdep) ) + ALLOCATE( champ_msk(imdep*jmdep) ) + +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,latmsk_ini) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,latmsk_ini) +#endif +c + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,varid,champ_msk) +#else + ierr = NF_GET_VAR_REAL(ncid,varid,champ_msk) +#endif +c + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF +c + title='RELIEF' + + CALL conf_dat2d(title,imdep, jmdep, lonmsk_ini, latmsk_ini, + . dlon_msk, dlat_msk, champ_msk, interbar ) + + CALL mask_c_o(imdep, jmdep, dlon_msk, dlat_msk,champ_msk, + . iim, jjp1, rlonv, rlatu, champint) + CALL gr_int_dyn(champint, masque, iim, jjp1) + DO i = 1, iim + masque(i,1) = FLOAT(NINT(masque(i,1))) + masque(i,jjp1) = FLOAT(NINT(masque(i,jjp1))) + ENDDO + DO i = 1, iim + DO j = 1, jjp1 + mask(i,j) = champint(i,j) + ENDDO + ENDDO + CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, phy_nat0) + ierr = NF_CLOSE(ncid) +c +c +C Traitement de la rugosite +c + PRINT*, 'Traitement de la rugosite' + ierr = NF_OPEN('Rugos.nc', NF_NOWRITE, ncid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ierr = NF_INQ_VARID(ncid,'RUGOS',varid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_VARDIMID (ncid,varid,ndimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, 'dimension ', imdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ALLOCATE( dlon_ini(imdep) ) + ALLOCATE( dlon(imdep) ) + +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, 'dimension ', jmdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ALLOCATE( dlat_ini(jmdep) ) + ALLOCATE( dlat(jmdep) ) + +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, 'dimension ', lmdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF +c + ALLOCATE( champ(imdep*jmdep) ) + + DO 200 l = 1, lmdep + dimfirst(1) = 1 + dimfirst(2) = 1 + dimfirst(3) = l +c + dimlast(1) = imdep + dimlast(2) = jmdep + dimlast(3) = 1 +c + PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l) + print*,dimfirst,dimlast +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ) +#else + ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + title = 'Rugosite Amip ' +c + CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini, + . dlon, dlat, champ, interbar ) + + IF ( interbar ) THEN + DO j = 1, imdep * jmdep + champ(j) = LOG(champ(j)) + ENDDO + + IF( l.EQ.1 ) THEN + WRITE(6,*) '-------------------------------------------------', + ,'------------------------' + WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', + , ' pour la rugosite $$$ ' + WRITE(6,*) '-------------------------------------------------', + ,'------------------------' + ENDIF + CALL inter_barxy ( imdep,jmdep -1,dlon,dlat,champ , + , iim,jjm,rlonu,rlatv, jjp1,champint ) + DO j=1,jjp1 + DO i=1,iim + champint(i,j)=EXP(champint(i,j)) + ENDDO + ENDDO + + DO j = 1, jjp1 + DO i = 1, iim + IF(NINT(mask(i,j)).NE.1) THEN + champint( i,j ) = 0.001 + ENDIF + ENDDO + ENDDO + ELSE + CALL rugosite(imdep, jmdep, dlon, dlat, champ, + . iim, jjp1, rlonv, rlatu, champint, mask) + ENDIF + DO j = 1,jjp1 + DO i = 1, iim + champtime (i,j,l) = champint(i,j) + ENDDO + ENDDO +200 CONTINUE +c + DO l = 1, lmdep + timeyear(l) = timecoord(l) + ENDDO + + PRINT 222, timeyear +222 FORMAT(2x,' Time year ',10f6.1) +c + + PRINT*, 'Interpolation temporelle dans l annee' + + DO j = 1, jjp1 + DO i = 1, iim + DO l = 1, lmdep + ax(l) = timeyear(l) + ay(l) = champtime (i,j,l) + ENDDO + CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder) + DO k = 1, 360 + time = FLOAT(k-1) + CALL SPLINT(ax,ay,yder,lmdep,time,by) + champan(i,j,k) = by + ENDDO + ENDDO + ENDDO + DO k = 1, 360 + DO j = 1, jjp1 + champan(iip1,j,k) = champan(1,j,k) + ENDDO + IF ( k.EQ.10 ) THEN + DO j = 1, jjp1 + CALL minmax( iip1,champan(1,j,10),chmin,chmax ) + PRINT *,' Rugosite au temps 10 ', chmin,chmax,j + ENDDO + ENDIF + ENDDO +c + DO k = 1, 360 + CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k), phy_rug(1,k)) + ENDDO +c + ierr = NF_CLOSE(ncid) + + DEALLOCATE( dlon ) + DEALLOCATE( dlon_ini ) + DEALLOCATE( dlat ) + DEALLOCATE( dlat_ini ) + DEALLOCATE( champ ) +c +c +C Traitement de la glace oceanique +c + PRINT*, 'Traitement de la glace oceanique' + ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ierr = NF_INQ_VARID(ncid,'SEA_ICE',varid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_VARDIMID (ncid,varid,ndimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, 'dimension ', imdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ALLOCATE ( dlon_ini(imdep) ) + ALLOCATE ( dlon(imdep) ) + +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, jmdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ALLOCATE ( dlat_ini(jmdep) ) + ALLOCATE ( dlat(jmdep) ) + +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, lmdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF +c + ALLOCATE ( champ(imdep*jmdep) ) + + DO l = 1, lmdep + dimfirst(1) = 1 + dimfirst(2) = 1 + dimfirst(3) = l +c + dimlast(1) = imdep + dimlast(2) = jmdep + dimlast(3) = 1 +c + PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l) +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ) +#else + ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + title = 'Sea-ice Amip ' +c + CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini, + . dlon, dlat, champ, interbar ) +c + IF( oldice ) THEN + CALL sea_ice(imdep, jmdep, dlon, dlat, champ, + . iim, jjp1, rlonv, rlatu, champint ) + ELSEIF ( interbar ) THEN + IF( l.EQ.1 ) THEN + WRITE(6,*) '-------------------------------------------------', + ,'------------------------' + WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', + , ' pour Sea-ice Amip $$$ ' + WRITE(6,*) '-------------------------------------------------', + ,'------------------------' + ENDIF + + CALL inter_barxy ( imdep,jmdep -1,dlon, dlat , + , champ, iim, jjm, rlonu, rlatv, jjp1, champint ) + ELSE + CALL sea_ice(imdep, jmdep, dlon, dlat, champ, + . iim, jjp1, rlonv, rlatu, champint ) + ENDIF + DO j = 1,jjp1 + DO i = 1, iim + champtime (i,j,l) = champint(i,j) + ENDDO + ENDDO + ENDDO +c + DO l = 1, lmdep + timeyear(l) = timecoord(l) + ENDDO + PRINT 222, timeyear +c + PRINT*, 'Interpolation temporelle' + DO j = 1, jjp1 + DO i = 1, iim + DO l = 1, lmdep + ax(l) = timeyear(l) + ay(l) = champtime (i,j,l) + ENDDO + CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder) + DO k = 1, 360 + time = FLOAT(k-1) + CALL SPLINT(ax,ay,yder,lmdep,time,by) + champan(i,j,k) = by + ENDDO + ENDDO + ENDDO + DO k = 1, 360 + DO j = 1, jjp1 + champan(iip1, j, k) = champan(1, j, k) + ENDDO + IF ( k.EQ.10 ) THEN + DO j = 1, jjp1 + CALL minmax( iip1,champan(1,j,10),chmin,chmax ) + PRINT *,' Sea ice au temps 10 ', chmin,chmax,j + ENDDO + ENDIF + ENDDO +c + + DO k = 1, 360 + CALL gr_dyn_fi(1, iip1, jjp1, klon, + . champan(1,1,k), phy_ice(1,k)) + DO i = 1, klon + phy_nat(i,k) = phy_nat0(i) + IF ( (phy_ice(i,k) - 0.5).GE.1.e-5 ) THEN + IF (NINT(phy_nat0(i)).EQ.0) THEN + phy_nat(i,k) = 3.0 + ELSE + phy_nat(i,k) = 2.0 + ENDIF + ENDIF + IF( NINT(phy_nat(i,k)).EQ.0 ) THEN + IF ( phy_rug(i,k).NE.0.001 ) phy_rug(i,k) = 0.001 + ENDIF + ENDDO + + ENDDO +c + + ierr = NF_CLOSE(ncid) +c + DEALLOCATE( dlon ) + DEALLOCATE( dlon_ini ) + DEALLOCATE( dlat ) + DEALLOCATE( dlat_ini ) + DEALLOCATE( champ ) + +477 continue +c +C Traitement de la sst +c + PRINT*, 'Traitement de la sst' + ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ierr = NF_INQ_VARID(ncid,'SST',varid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_VARDIMID (ncid,varid,ndimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable SST ', namedim,'dimension ', imdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ALLOCATE( dlon_ini(imdep) ) + ALLOCATE( dlon(imdep) ) + +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini) +#endif + + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable SST ', namedim, 'dimension ', jmdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ALLOCATE( dlat_ini(jmdep) ) + ALLOCATE( dlat(jmdep) ) + +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, 'dimension ', lmdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ALLOCATE( champ(imdep*jmdep) ) + IF( extrap ) THEN + ALLOCATE ( work(imdep,jmdep) ) + ENDIF +c + DO l = 1, lmdep + dimfirst(1) = 1 + dimfirst(2) = 1 + dimfirst(3) = l +c + dimlast(1) = imdep + dimlast(2) = jmdep + dimlast(3) = 1 +c + PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l) +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ) +#else + ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + title='Sst Amip' +c + CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini, + . dlon, dlat, champ, interbar ) + IF ( extrap ) THEN + CALL extrapol(champ, imdep, jmdep, 999999.,.TRUE.,.TRUE.,2,work) + ENDIF +c + + IF ( interbar ) THEN + IF( l.EQ.1 ) THEN + WRITE(6,*) '-------------------------------------------------', + ,'------------------------' + WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', + , ' pour la Sst Amip $$$ ' + WRITE(6,*) '-------------------------------------------------', + ,'------------------------' + ENDIF + CALL inter_barxy ( imdep,jmdep -1,dlon, dlat , + , champ, iim, jjm, rlonu, rlatv, jjp1, champint ) + ELSE + CALL grille_m(imdep, jmdep, dlon, dlat, champ, + . iim, jjp1, rlonv, rlatu, champint ) + ENDIF + + DO j = 1,jjp1 + DO i = 1, iim + champtime (i,j,l) = champint(i,j) + ENDDO + ENDDO + ENDDO +c + DO l = 1, lmdep + timeyear(l) = timecoord(l) + ENDDO + print 222, timeyear +c +C interpolation temporelle + DO j = 1, jjp1 + DO i = 1, iim + DO l = 1, lmdep + ax(l) = timeyear(l) + ay(l) = champtime (i,j,l) + ENDDO + CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder) + DO k = 1, 360 + time = FLOAT(k-1) + CALL SPLINT(ax,ay,yder,lmdep,time,by) + champan(i,j,k) = by + ENDDO + ENDDO + ENDDO + DO k = 1, 360 + DO j = 1, jjp1 + champan(iip1,j,k) = champan(1,j,k) + ENDDO + IF ( k.EQ.10 ) THEN + DO j = 1, jjp1 + CALL minmax( iip1,champan(1,j,10),chmin,chmax ) + PRINT *,' SST au temps 10 ', chmin,chmax,j + ENDDO + ENDIF + ENDDO +c + DO k = 1, 360 + CALL gr_dyn_fi(1, iip1, jjp1, klon, + . champan(1,1,k), phy_sst(1,k)) + ENDDO +c + ierr = NF_CLOSE(ncid) +c + DEALLOCATE( dlon ) + DEALLOCATE( dlon_ini ) + DEALLOCATE( dlat ) + DEALLOCATE( dlat_ini ) + DEALLOCATE( champ ) +c +C Traitement de l'albedo +c + PRINT*, 'Traitement de l albedo' + ierr = NF_OPEN('Albedo.nc', NF_NOWRITE, ncid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_VARID(ncid,'ALBEDO',varid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_VARDIMID (ncid,varid,ndimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, 'dimension ', imdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ALLOCATE ( dlon_ini(imdep) ) + ALLOCATE ( dlon(imdep) ) + +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, 'dimension ', jmdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + ALLOCATE ( dlat_ini(jmdep) ) + ALLOCATE ( dlat(jmdep) ) + +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + print*,'variable ', namedim, 'dimension ', lmdep + ierr = NF_INQ_VARID(ncid,namedim,dimid) + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord) +#else + ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF +c + ALLOCATE ( champ(imdep*jmdep) ) + + DO l = 1, lmdep + dimfirst(1) = 1 + dimfirst(2) = 1 + dimfirst(3) = l +c + dimlast(1) = imdep + dimlast(2) = jmdep + dimlast(3) = 1 +c + PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l) +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ) +#else + ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ) +#endif + if (ierr.ne.0) then + print *, NF_STRERROR(ierr) + STOP + ENDIF + + title='Albedo Amip' +c + CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini, + . dlon, dlat, champ, interbar ) +c +c + IF ( interbar ) THEN + IF( l.EQ.1 ) THEN + WRITE(6,*) '-------------------------------------------------', + ,'------------------------' + WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', + , ' pour l Albedo Amip $$$ ' + WRITE(6,*) '-------------------------------------------------', + ,'------------------------' + ENDIF + + CALL inter_barxy ( imdep,jmdep -1,dlon, dlat , + , champ, iim, jjm, rlonu, rlatv, jjp1, champint ) + ELSE + CALL grille_m(imdep, jmdep, dlon, dlat, champ, + . iim, jjp1, rlonv, rlatu, champint ) + ENDIF +c + DO j = 1,jjp1 + DO i = 1, iim + champtime (i, j, l) = champint(i, j) + ENDDO + ENDDO + ENDDO +c + DO l = 1, lmdep + timeyear(l) = timecoord(l) + ENDDO + print 222, timeyear +c +C interpolation temporelle + DO j = 1, jjp1 + DO i = 1, iim + DO l = 1, lmdep + ax(l) = timeyear(l) + ay(l) = champtime (i, j, l) + ENDDO + CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder) + DO k = 1, 360 + time = FLOAT(k-1) + CALL SPLINT(ax,ay,yder,lmdep,time,by) + champan(i,j,k) = by + ENDDO + ENDDO + ENDDO + DO k = 1, 360 + DO j = 1, jjp1 + champan(iip1, j, k) = champan(1, j, k) + ENDDO + IF ( k.EQ.10 ) THEN + DO j = 1, jjp1 + CALL minmax( iip1,champan(1,j,10),chmin,chmax ) + PRINT *,' Albedo au temps 10 ', chmin,chmax,j + ENDDO + ENDIF + ENDDO +c + DO k = 1, 360 + CALL gr_dyn_fi(1, iip1, jjp1, klon, + . champan(1,1,k), phy_alb(1,k)) + ENDDO +c + ierr = NF_CLOSE(ncid) +c +c + DO k = 1, 360 + DO i = 1, klon + phy_bil(i,k) = 0.0 + ENDDO + ENDDO +c + PRINT*, 'Ecriture du fichier limit' +c + ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid) +c + ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30, + . "Fichier conditions aux limites") + ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim) + ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim) +c + dims(1) = ndim + dims(2) = ntim +c +ccc ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim) + ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim) + ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17, + . "Jour dans l annee") +ccc ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT) + ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT) + ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23, + . "Nature du sol (0,1,2,3)") +ccc ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST) + ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST) + ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35, + . "Temperature superficielle de la mer") +ccc ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS) + ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS) + ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32, + . "Reference flux de chaleur au sol") +ccc ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB) + ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB) + ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19, + . "Albedo a la surface") +ccc ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG) + ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG) + ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8, + . "Rugosite") +c + ierr = NF_ENDDEF(nid) +c + DO k = 1, 360 +c + debut(1) = 1 + debut(2) = k + epais(1) = klon + epais(2) = 1 +c +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k)) +#else + ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k)) + ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k)) +#endif +c + ENDDO +c + ierr = NF_CLOSE(nid) +c + STOP + END Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/sort.F =================================================================== --- LMDZ.3.3/branches/rel-LF/libf/dyn3d/sort.F (revision 265) +++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/sort.F (revision 265) @@ -0,0 +1,35 @@ +C +C $Header$ +C + SUBROUTINE sort(n,d) +c +c P.Le Van +c +c... cette routine met le tableau d dans l'ordre croissant .... +cc ( pour avoir l'ordre decroissant,il suffit de remplacer l'instruc +c tion situee + bas IF(d(j).LE.p) THEN par +c IF(d(j).GE.p) THEN +c + + INTEGER n + REAL d(1) , p + INTEGER i,j,k + + DO i=1,n-1 + k=i + p=d(i) + DO j=i+1,n + IF(d(j).LE.p) THEN + k=j + p=d(j) + ENDIF + ENDDO + + IF(k.ne.i) THEN + d(k)=d(i) + d(i)=p + ENDIF + ENDDO + + RETURN + END Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/conemav.F =================================================================== --- LMDZ.3.3/branches/rel-LF/libf/phylmd/conemav.F (revision 265) +++ LMDZ.3.3/branches/rel-LF/libf/phylmd/conemav.F (revision 265) @@ -0,0 +1,147 @@ + SUBROUTINE conemav (dtime,paprs,pplay,t,q,u,v,tra,ntra, + . work1,work2,d_t,d_q,d_u,d_v,d_tra, + . rain, snow, kbas, ktop, + . upwd,dnwd,dnwdbis,Ma,cape,tvp,iflag, + . pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr) + +c + IMPLICIT none +c====================================================================== +c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818 +c Objet: schema de convection de Emanuel (1991) interface +c====================================================================== +c Arguments: +c dtime--input-R-pas d'integration (s) +c s-------input-R-la valeur "s" pour chaque couche +c sigs----input-R-la valeur "sigma" de chaque couche +c sig-----input-R-la valeur de "sigma" pour chaque niveau +c psolpa--input-R-la pression au sol (en Pa) +C pskapa--input-R-exponentiel kappa de psolpa +c h-------input-R-enthalpie potentielle (Cp*T/P**kappa) +c q-------input-R-vapeur d'eau (en kg/kg) +c +c work*: input et output: deux variables de travail, +c on peut les mettre a 0 au debut +c ALE-----input-R-energie disponible pour soulevement +c +C d_h-----output-R-increment de l'enthalpie potentielle (h) +c d_q-----output-R-increment de la vapeur d'eau +c rain----output-R-la pluie (mm/s) +c snow----output-R-la neige (mm/s) +c upwd----output-R-saturated updraft mass flux (kg/m**2/s) +c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s) +c dnwd0---output-R-unsaturated downdraft mass flux (kg/m**2/s) +c Cape----output-R-CAPE (J/kg) +c Tvp-----output-R-Temperature virtuelle d'une parcelle soulevee +c adiabatiquement a partir du niveau 1 (K) +c deltapb-output-R-distance entre LCL et base de la colonne (<0 ; Pa) +c Ice_flag-input-L-TRUE->prise en compte de la thermodynamique de la glace +c====================================================================== +c +#include "dimensions.h" +#include "dimphy.h" +c + integer NTRAC + PARAMETER (NTRAC=nqmx-2) +c + REAL dtime, paprs(klon,klev+1),pplay(klon,klev) + REAL t(klon,klev),q(klon,klev),u(klon,klev),v(klon,klev) + REAL tra(klon,klev,ntrac) + INTEGER ntra + REAL work1(klon,klev),work2(klon,klev) +c + REAL d_t(klon,klev),d_q(klon,klev),d_u(klon,klev),d_v(klon,klev) + REAL d_tra(klon,klev,ntrac) + REAL rain(klon),snow(klon) +c + INTEGER kbas(klon),ktop(klon) + REAL em_ph(klon,klev+1),em_p(klon,klev) + REAL upwd(klon,klev),dnwd(klon,klev),dnwdbis(klon,klev) + REAL Ma(klon,klev),cape(klon),tvp(klon,klev) + INTEGER iflag(klon) + REAL rflag(klon) + REAL pbase(klon),bbase(klon) + REAL dtvpdt1(klon,klev),dtvpdq1(klon,klev) + REAL dplcldt(klon),dplcldr(klon) +c + REAL zx_t,zdelta,zx_qs,zcor +c + INTEGER noff, minorig + INTEGER i,k,itra + REAL qs(klon,klev) + REAL cbmf(klon) + SAVE cbmf + INTEGER ifrst + SAVE ifrst + DATA ifrst /0/ +#include "YOMCST.h" +#include "YOETHF.h" +#include "FCTTRE.h" +c +c + IF (ifrst .EQ. 0) THEN + ifrst = 1 + DO i = 1, klon + cbmf(i) = 0. + ENDDO + ENDIF + + DO k = 1, klev+1 + DO i=1,klon + em_ph(i,k) = paprs(i,k) / 100.0 + ENDDO + ENDDO +c + DO k = 1, klev + DO i=1,klon + em_p(i,k) = pplay(i,k) / 100.0 + ENDDO + ENDDO + +c + DO k = 1, klev + DO i = 1, klon + zx_t = t(i,k) + zdelta=MAX(0.,SIGN(1.,rtt-zx_t)) + zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0) + zcor=1./(1.-retv*zx_qs) + qs(i,k)=zx_qs*zcor + ENDDO + ENDDO +c + noff = 2 + minorig = 2 + CALL convect1(klon,klev,klev+1,noff,minorig,t,q,qs,u,v, + $ em_p,em_ph,iflag, + $ d_t,d_q,d_u,d_v,rain,cbmf,dtime,Ma) +c + DO i = 1,klon + rain(i) = rain(i)/86400. + rflag(i)=iflag(i) + ENDDO +c call dump2d(iim,jjm-1,rflag(2:klon-1),'FLAG CONVECTION ') +c if (klon.eq.1) then +c print*,'IFLAG ',iflag +c else +c write(*,'(96i1)') (iflag(i),i=2,klon-1) +c endif + DO k = 1, klev + DO i = 1, klon + d_t(i,k) = dtime*d_t(i,k) + d_q(i,k) = dtime*d_q(i,k) + d_u(i,k) = dtime*d_u(i,k) + d_v(i,k) = dtime*d_v(i,k) + ENDDO + DO itra = 1,ntra + DO i = 1, klon + d_tra(i,k,itra) = 0. + ENDDO + ENDDO + ENDDO + +c +c +c + RETURN + END + Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/convect1.F =================================================================== --- LMDZ.3.3/branches/rel-LF/libf/phylmd/convect1.F (revision 265) +++ LMDZ.3.3/branches/rel-LF/libf/phylmd/convect1.F (revision 265) @@ -0,0 +1,645 @@ + subroutine convect1(len,nd,ndp1,noff,minorig, + & t,q,qs,u,v, + & p,ph,iflag,ft, + & fq,fu,fv,precip,cbmf,delt,Ma) +C.............................START PROLOGUE............................ +C +C SCCS IDENTIFICATION: @(#)convect1.f 1.1 04/21/00 +C 19:40:52 /h/cm/library/nogaps4/src/sub/fcst/convect1.f_v +C +C CONFIGURATION IDENTIFICATION: None +C +C MODULE NAME: convect1 +C +C DESCRIPTION: +C +C convect1 The Emanuel Cumulus Convection Scheme +C +C CONTRACT NUMBER AND TITLE: None +C +C REFERENCES: Programmers K. Emanuel (MIT), Timothy F. Hogan, M. Peng (NRL) +C +C CLASSIFICATION: Unclassified +C +C RESTRICTIONS: None +C +C COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90 +C +C COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress" +C Fortran 90: -O vector3,scalar3,task1,aggress,overindex -ei -r 2 +C +C LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a +C +C USAGE: call convect1(len,nd,noff,minorig, +C & t,q,qs,u,v, +C & p,ph,iflag,ft, +C & fq,fu,fv,precip,cbmf,delt) +C +C PARAMETERS: +C Name Type Usage Description +C ---------- ---------- ------- ---------------------------- +C +C len Integer Input first (i) dimension +C nd Integer Input vertical (k) dimension +C ndp1 Integer Input nd + 1 +C noff Integer Input integer limit for convection (nd-noff) +C minorig Integer Input First level of convection +C t Real Input temperature +C q Real Input specific hum +C qs Real Input sat specific hum +C u Real Input u-wind +C v Real Input v-wind +C p Real Input full level pressure +C ph Real Input half level pressure +C iflag Integer Output iflag on latitude strip +C ft Real Output temp tend +C fq Real Output spec hum tend +C fu Real Output u-wind tend +C fv Real Output v-wind tend +C cbmf Real In/Out cumulus mass flux +C delt Real Input time step +C iflag Integer Output integer flag for Emanuel conditions +C +C COMMON BLOCKS: +C Block Name Type Usage Notes +C -------- -------- ---- ------ ------------------------ +C +C FILES: None +C +C DATA BASES: None +C +C NON-FILE INPUT/OUTPUT: None +C +C ERROR CONDITIONS: None +C +C ADDITIONAL COMMENTS: None +C +C.................MAINTENANCE SECTION................................ +C +C MODULES CALLED: +C Name Description +C convect2 Emanuel cumulus convection tendency calculations +C ------- ---------------------- +C LOCAL VARIABLES AND +C STRUCTURES: +C Name Type Description +C ------- ------ ----------- +C See Comments Below +C +C i Integer loop index +C k Integer loop index +c +C METHOD: +C +C See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation of a +C convective scheme for use in climate models. +C +C FILES: None +C +C INCLUDE FILES: None +C +C MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak +C +C..............................END PROLOGUE............................. +c +c + implicit none +c +#include "dimensions.h" +#include "dimphy.h" +c + integer len + integer nd + integer ndp1 + integer noff + real t(len,nd) + real q(len,nd) + real qs(len,nd) + real u(len,nd) + real v(len,nd) + real p(len,nd) + real ph(len,ndp1) + integer iflag(len) + real ft(len,nd) + real fq(len,nd) + real fu(len,nd) + real fv(len,nd) + real precip(len) + real cbmf(len) + real Ma(len,nd) + integer minorig + real delt,cpd,cpv,cl,rv,rd,lv0,g + real sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp + real alpha,entp,coeffs,coeffr,omtrain,cu +c +!------------------------------------------------------------------- +! --- ARGUMENTS +!------------------------------------------------------------------- +! --- On input: +! +! t: Array of absolute temperature (K) of dimension ND, with first +! index corresponding to lowest model level. Note that this array +! will be altered by the subroutine if dry convective adjustment +! occurs and if IPBL is not equal to 0. +! +! q: Array of specific humidity (gm/gm) of dimension ND, with first +! index corresponding to lowest model level. Must be defined +! at same grid levels as T. Note that this array will be altered +! if dry convective adjustment occurs and if IPBL is not equal to 0. +! +! qs: Array of saturation specific humidity of dimension ND, with first +! index corresponding to lowest model level. Must be defined +! at same grid levels as T. Note that this array will be altered +! if dry convective adjustment occurs and if IPBL is not equal to 0. +! +! u: Array of zonal wind velocity (m/s) of dimension ND, witth first +! index corresponding with the lowest model level. Defined at +! same levels as T. Note that this array will be altered if +! dry convective adjustment occurs and if IPBL is not equal to 0. +! +! v: Same as u but for meridional velocity. +! +! tra: Array of passive tracer mixing ratio, of dimensions (ND,NTRA), +! where NTRA is the number of different tracers. If no +! convective tracer transport is needed, define a dummy +! input array of dimension (ND,1). Tracers are defined at +! same vertical levels as T. Note that this array will be altered +! if dry convective adjustment occurs and if IPBL is not equal to 0. +! +! p: Array of pressure (mb) of dimension ND, with first +! index corresponding to lowest model level. Must be defined +! at same grid levels as T. +! +! ph: Array of pressure (mb) of dimension ND+1, with first index +! corresponding to lowest level. These pressures are defined at +! levels intermediate between those of P, T, Q and QS. The first +! value of PH should be greater than (i.e. at a lower level than) +! the first value of the array P. +! +! nl: The maximum number of levels to which convection can penetrate, plus 1. +! NL MUST be less than or equal to ND-1. +! +! delt: The model time step (sec) between calls to CONVECT +! +!---------------------------------------------------------------------------- +! --- On Output: +! +! iflag: An output integer whose value denotes the following: +! VALUE INTERPRETATION +! ----- -------------- +! 0 Moist convection occurs. +! 1 Moist convection occurs, but a CFL condition +! on the subsidence warming is violated. This +! does not cause the scheme to terminate. +! 2 Moist convection, but no precip because ep(inb) lt 0.0001 +! 3 No moist convection because new cbmf is 0 and old cbmf is 0. +! 4 No moist convection; atmosphere is not +! unstable +! 6 No moist convection because ihmin le minorig. +! 7 No moist convection because unreasonable +! parcel level temperature or specific humidity. +! 8 No moist convection: lifted condensation +! level is above the 200 mb level. +! 9 No moist convection: cloud base is higher +! then the level NL-1. +! +! ft: Array of temperature tendency (K/s) of dimension ND, defined at same +! grid levels as T, Q, QS and P. +! +! fq: Array of specific humidity tendencies ((gm/gm)/s) of dimension ND, +! defined at same grid levels as T, Q, QS and P. +! +! fu: Array of forcing of zonal velocity (m/s^2) of dimension ND, +! defined at same grid levels as T. +! +! fv: Same as FU, but for forcing of meridional velocity. +! +! ftra: Array of forcing of tracer content, in tracer mixing ratio per +! second, defined at same levels as T. Dimensioned (ND,NTRA). +! +! precip: Scalar convective precipitation rate (mm/day). +! +! wd: A convective downdraft velocity scale. For use in surface +! flux parameterizations. See convect.ps file for details. +! +! tprime: A convective downdraft temperature perturbation scale (K). +! For use in surface flux parameterizations. See convect.ps +! file for details. +! +! qprime: A convective downdraft specific humidity +! perturbation scale (gm/gm). +! For use in surface flux parameterizations. See convect.ps +! file for details. +! +! cbmf: The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST +! BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT +! ITS NEXT CALL. That is, the value of CBMF must be "remembered" +! by the calling program between calls to CONVECT. +! +! det: Array of detrainment mass flux of dimension ND. +! +!------------------------------------------------------------------- +c +c Local arrays +c + integer nl + integer nlp + integer nlm + integer i,k,n + real delti + real rowl + real clmcpv + real clmcpd + real cpdmcp + real cpvmcpd + real eps + real epsi + real epsim1 + real ginv + real hrd + real prccon1 + integer icbmax + real lv(klon,klev) + real cpn(klon,klev) + real cpx(klon,klev) + real tv(klon,klev) + real gz(klon,klev) + real hm(klon,klev) + real h(klon,klev) + real work(klon) + integer ihmin(klon) + integer nk(klon) + real rh(klon) + real chi(klon) + real plcl(klon) + integer icb(klon) + real tnk(klon) + real qnk(klon) + real gznk(klon) + real pnk(klon) + real qsnk(klon) + real ticb(klon) + real gzicb(klon) + real tp(klon,klev) + real tvp(klon,klev) + real clw(klon,klev) +c + real ah0(klon),cpp(klon) + real tg,qg,s,alv,tc,ahg,denom,es,rg +c + integer ncum + integer idcum(klon) +c + cpd=1005.7 + cpv=1870.0 + cl=4190.0 + rv=461.5 + rd=287.04 + lv0=2.501E6 + g=9.8 +C +C *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) *** +C *** TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- *** +C *** CONVERSION THRESHOLD IS ASSUMED TO BE ZERO *** +C *** (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY *** +C *** BETWEEN 0 C AND TLCRIT) *** +C *** ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT *** +C *** FORMULATION *** +C *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** +C *** SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** +C *** OF CLOUD *** +C *** OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN *** +C *** OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW *** +C *** COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** +C *** OF RAIN *** +C *** COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** +C *** OF SNOW *** +C *** CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM *** +C *** TRANSPORT *** +C *** DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION *** +C *** A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC *** +C *** ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF *** +C *** APPROACH TO QUASI-EQUILIBRIUM *** +C *** (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) *** +C *** (DAMP MUST BE LESS THAN 1) *** +c + sigs=0.12 + sigd=0.05 + elcrit=0.0011 + tlcrit=-55.0 + omtsnow=5.5 + dtmax=0.9 + damp=0.1 + alpha=0.2 + entp=1.5 + coeffs=0.8 + coeffr=1.0 + omtrain=50.0 +c + cu=0.70 + damp=0.1 +c +c +c Define nl, nlp, nlm, and delti +c + nl=nd-noff + nlp=nl+1 + nlm=nl-1 + delti=1.0/delt +! +!------------------------------------------------------------------- +! --- SET CONSTANTS +!------------------------------------------------------------------- +! + rowl=1000.0 + clmcpv=cl-cpv + clmcpd=cl-cpd + cpdmcp=cpd-cpv + cpvmcpd=cpv-cpd + eps=rd/rv + epsi=1.0/eps + epsim1=epsi-1.0 + ginv=1.0/g + hrd=0.5*rd + prccon1=86400.0*1000.0/(rowl*g) +! +! dtmax is the maximum negative temperature perturbation. +! +!===================================================================== +! --- INITIALIZE OUTPUT ARRAYS AND PARAMETERS +!===================================================================== +! + do 20 k=1,nd + do 10 i=1,len + ft(i,k)=0.0 + fq(i,k)=0.0 + fu(i,k)=0.0 + fv(i,k)=0.0 + tvp(i,k)=0.0 + tp(i,k)=0.0 + clw(i,k)=0.0 + gz(i,k) = 0. + 10 continue + 20 continue + do 60 i=1,len + precip(i)=0.0 + iflag(i)=0 + 60 continue +c +!===================================================================== +! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY +!===================================================================== + do 110 k=1,nl+1 + do 100 i=1,len + lv(i,k)= lv0-clmcpv*(t(i,k)-273.15) + cpn(i,k)=cpd*(1.0-q(i,k))+cpv*q(i,k) + cpx(i,k)=cpd*(1.0-q(i,k))+cl*q(i,k) + tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1) + 100 continue + 110 continue +c +c gz = phi at the full levels (same as p). +c + do 120 i=1,len + gz(i,1)=0.0 + 120 continue + do 140 k=2,nlp + do 130 i=1,len + gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k)) + & *(p(i,k-1)-p(i,k))/ph(i,k) + 130 continue + 140 continue +c +c h = phi + cpT (dry static energy). +c hm = phi + cp(T-Tbase)+Lq +c + do 170 k=1,nlp + do 160 i=1,len + h(i,k)=gz(i,k)+cpn(i,k)*t(i,k) + hm(i,k)=gz(i,k)+cpx(i,k)*(t(i,k)-t(i,1))+lv(i,k)*q(i,k) + 160 continue + 170 continue +c +!------------------------------------------------------------------- +! --- Find level of minimum moist static energy +! --- If level of minimum moist static energy coincides with +! --- or is lower than minimum allowable parcel origin level, +! --- set iflag to 6. +!------------------------------------------------------------------- + do 180 i=1,len + work(i)=1.0e12 + ihmin(i)=nl + 180 continue + do 200 k=2,nlp + do 190 i=1,len + if((hm(i,k).lt.work(i)).and. + & (hm(i,k).lt.hm(i,k-1)))then + work(i)=hm(i,k) + ihmin(i)=k + endif + 190 continue + 200 continue + do 210 i=1,len + ihmin(i)=min(ihmin(i),nlm) + if(ihmin(i).le.minorig)then + iflag(i)=6 + endif + 210 continue +c +!------------------------------------------------------------------- +! --- Find that model level below the level of minimum moist static +! --- energy that has the maximum value of moist static energy +!------------------------------------------------------------------- + + do 220 i=1,len + work(i)=hm(i,minorig) + nk(i)=minorig + 220 continue + do 240 k=minorig+1,nl + do 230 i=1,len + if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then + work(i)=hm(i,k) + nk(i)=k + endif + 230 continue + 240 continue +!------------------------------------------------------------------- +! --- Check whether parcel level temperature and specific humidity +! --- are reasonable +!------------------------------------------------------------------- + do 250 i=1,len + if(((t(i,nk(i)).lt.250.0).or. + & (q(i,nk(i)).le.0.0).or. + & (p(i,ihmin(i)).lt.400.0)).and. + & (iflag(i).eq.0))iflag(i)=7 + 250 continue +!------------------------------------------------------------------- +! --- Calculate lifted condensation level of air at parcel origin level +! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980) +!------------------------------------------------------------------- + do 260 i=1,len + tnk(i)=t(i,nk(i)) + qnk(i)=q(i,nk(i)) + gznk(i)=gz(i,nk(i)) + pnk(i)=p(i,nk(i)) + qsnk(i)=qs(i,nk(i)) +c + rh(i)=qnk(i)/qsnk(i) + rh(i)=min(1.0,rh(i)) + chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i)) + plcl(i)=pnk(i)*(rh(i)**chi(i)) + if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0)) + & .and.(iflag(i).eq.0))iflag(i)=8 + 260 continue +!------------------------------------------------------------------- +! --- Calculate first level above lcl (=icb) +!------------------------------------------------------------------- + do 270 i=1,len + icb(i)=nlm + 270 continue +c + do 290 k=minorig,nl + do 280 i=1,len + if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i))) + & icb(i)=min(icb(i),k) + 280 continue + 290 continue +c + do 300 i=1,len + if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9 + 300 continue +c +c Compute icbmax. +c + icbmax=2 + do 310 i=1,len + icbmax=max(icbmax,icb(i)) + 310 continue +! +!------------------------------------------------------------------- +! --- Calculates the lifted parcel virtual temperature at nk, +! --- the actual temperature, and the adiabatic +! --- liquid water content. The procedure is to solve the equation. +! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. +!------------------------------------------------------------------- +! + do 320 i=1,len + tnk(i)=t(i,nk(i)) + qnk(i)=q(i,nk(i)) + gznk(i)=gz(i,nk(i)) + ticb(i)=t(i,icb(i)) + gzicb(i)=gz(i,icb(i)) + 320 continue +c +c *** Calculate certain parcel quantities, including static energy *** +c + do 330 i=1,len + ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + & +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i) + cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv + 330 continue +c +c *** Calculate lifted parcel quantities below cloud base *** +c + do 350 k=minorig,icbmax-1 + do 340 i=1,len + tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))/cpp(i) + tvp(i,k)=tp(i,k)*(1.+qnk(i)*epsi) + 340 continue + 350 continue +c +c *** Find lifted parcel quantities above cloud base *** +c + do 360 i=1,len + tg=ticb(i) + qg=qs(i,icb(i)) + alv=lv0-clmcpv*(ticb(i)-273.15) +c +c First iteration. +c + s=cpd+alv*alv*qg/(rv*ticb(i)*ticb(i)) + s=1./s + ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) + tg=tg+s*(ah0(i)-ahg) + tg=max(tg,35.0) + tc=tg-273.15 + denom=243.5+tc + if(tc.ge.0.0)then + es=6.112*exp(17.67*tc/denom) + else + es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) + endif + qg=eps*es/(p(i,icb(i))-es*(1.-eps)) +c +c Second iteration. +c + s=cpd+alv*alv*qg/(rv*ticb(i)*ticb(i)) + s=1./s + ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) + tg=tg+s*(ah0(i)-ahg) + tg=max(tg,35.0) + tc=tg-273.15 + denom=243.5+tc + if(tc.ge.0.0)then + es=6.112*exp(17.67*tc/denom) + else + es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) + end if + qg=eps*es/(p(i,icb(i))-es*(1.-eps)) +c + alv=lv0-clmcpv*(ticb(i)-273.15) + tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i) + & -gz(i,icb(i))-alv*qg)/cpd + clw(i,icb(i))=qnk(i)-qg + clw(i,icb(i))=max(0.0,clw(i,icb(i))) + rg=qg/(1.-qnk(i)) + tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi) + 360 continue +c + do 380 k=minorig,icbmax + do 370 i=1,len + tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i) + 370 continue + 380 continue +c +!------------------------------------------------------------------- +! --- Test for instability. +! --- If there was no convection at last time step and parcel +! --- is stable at icb, then set iflag to 4. +!------------------------------------------------------------------- + + do 390 i=1,len + if((cbmf(i).eq.0.0) .and.(iflag(i).eq.0).and. + & (tvp(i,icb(i)).le.(tv(i,icb(i))-dtmax)))iflag(i)=4 + 390 continue + +!===================================================================== +! --- IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY +!===================================================================== +c + ncum=0 + do 400 i=1,len + if(iflag(i).eq.0)then + ncum=ncum+1 + idcum(ncum)=i + endif + 400 continue +c +c Call convect2, which compresses the points and computes the heating, +c moistening, velocity mixing, and precipiation. +c +c print*,'cpd avant convect2 ',cpd + if(ncum.gt.0)then + call convect2(ncum,idcum,len,nd,ndp1,nl,minorig, + & nk,icb, + & t,q,qs,u,v,gz,tv,tp,tvp,clw,h, + & lv,cpn,p,ph,ft,fq,fu,fv, + & tnk,qnk,gznk,plcl, + & precip,cbmf,iflag, + & delt,cpd,cpv,cl,rv,rd,lv0,g, + & sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp, + & alpha,entp,coeffs,coeffr,omtrain,cu,Ma) + endif +c + return + end Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/convect2.F =================================================================== --- LMDZ.3.3/branches/rel-LF/libf/phylmd/convect2.F (revision 265) +++ LMDZ.3.3/branches/rel-LF/libf/phylmd/convect2.F (revision 265) @@ -0,0 +1,1391 @@ + subroutine convect2(ncum,idcum,len,nd,ndp1,nl,minorig, + & nk1,icb1, + & t1,q1,qs1,u1,v1,gz1,tv1,tp1,tvp1,clw1,h1, + & lv1,cpn1,p1,ph1,ft1,fq1,fu1,fv1, + & tnk1,qnk1,gznk1,plcl1, + & precip1,cbmf1,iflag1, + & delt,cpd,cpv,cl,rv,rd,lv0,g, + & sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp, + & alpha,entp,coeffs,coeffr,omtrain,cu,Ma) +C.............................START PROLOGUE............................ +C +C SCCS IDENTIFICATION: @(#)convect2.f 1.2 05/18/00 +C 22:06:22 /h/cm/library/nogaps4/src/sub/fcst/convect2.f_v +C +C CONFIGURATION IDENTIFICATION: None +C +C MODULE NAME: convect2 +C +C DESCRIPTION: +C +C convect1 The Emanuel Cumulus Convection Scheme - compute tendencies +C +C CONTRACT NUMBER AND TITLE: None +C +C REFERENCES: Programmers K. Emanuel (MIT), Timothy F. Hogan, M. Peng (NRL) +C +C CLASSIFICATION: Unclassified +C +C RESTRICTIONS: None +C +C COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90 +C +C COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress" +C Fortran 90: -O vector3,scalar3,task1,aggress,overindex -ei -r 2 +C +C LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a +C +C USAGE: call convect2(ncum,idcum,len,nd,nl,minorig, +C & nk1,icb1, +C & t1,q1,qs1,u1,v1,gz1,tv1,tp1,tvp1,clw1,h1, +C & lv1,cpn1,p1,ph1,ft1,fq1,fu1,fv1, +C & tnk1,qnk1,gznk1,plcl1, +C & precip1,cbmf1,iflag1, +C & delt,cpd,cpv,cl,rv,rd,lv0,g, +C & sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp, +C & alpha,entp,coeffs,coeffr,omtrain,cu) +C +C PARAMETERS: +C Name Type Usage Description +C ---------- ---------- ------- ---------------------------- +C +C ncum Integer Input number of cumulus points +C idcum Integer Input index of cumulus point +C len Integer Input first dimension +C nd Integer Input total vertical dimension +C ndp1 Integer Input nd + 1 +C nl Integer Input vertical dimension for convection +C minorig Integer Input First level where convection is allow to begin +C nk1 Integer Input First level of convection +C ncb1 Integer Input Level of free convection +C t1 Real Input temperature +C q1 Real Input specific hum +C qs1 Real Input sat specific hum +C u1 Real Input u-wind +C v1 Real Input v-wind +C gz1 Real Inout geop +C tv1 Real Input virtual temp +C tp1 Real Input +C clw1 Real Inout cloud liquid water +C h1 Real Inout +C lv1 Real Inout +C cpn1 Real Inout +C p1 Real Input full level pressure +C ph1 Real Input half level pressure +C ft1 Real Output temp tend +C fq1 Real Output spec hum tend +C fu1 Real Output u-wind tend +C fv1 Real Output v-wind tend +C precip1 Real Output prec +C cbmf1 Real In/Out cumulus mass flux +C iflag1 Integer Output iflag on latitude strip +C delt Real Input time step +C cpd Integer Input See description below +C cpv Integer Input See description below +C cl Integer Input See description below +C rv Integer Input See description below +C rd Integer Input See description below +C lv0 Integer Input See description below +C g Integer Input See description below +C sigs Integer Input See description below +C sigd Integer Input See description below +C elcrit Integer Input See description below +C tlcrit Integer Input See description below +C omtsnow Integer Input See description below +C dtmax Integer Input See description below +C damp Integer Input See description below +C alpha Integer Input See description below +C ent Integer Input See description below +C coeffs Integer Input See description below +C coeffr Integer Input See description below +C omtrain Integer Input See description below +C cu Integer Input See description below +C +C COMMON BLOCKS: +C Block Name Type Usage Notes +C -------- -------- ---- ------ ------------------------ +C +C FILES: None +C +C DATA BASES: None +C +C NON-FILE INPUT/OUTPUT: None +C +C ERROR CONDITIONS: None +C +C ADDITIONAL COMMENTS: None +C +C.................MAINTENANCE SECTION................................ +C +C MODULES CALLED: +C Name Description +C zilch Zero out an array +C ------- ---------------------- +C LOCAL VARIABLES AND +C STRUCTURES: +C Name Type Description +C ------- ------ ----------- +C See Comments Below +C +C i Integer loop index +C k Integer loop index +c +C METHOD: +C +C See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation of a +C convective scheme for use in climate models. +C +C FILES: None +C +C INCLUDE FILES: None +C +C MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak +C +C..............................END PROLOGUE............................. +c +c + implicit none +c +#include "dimensions.h" +#include "dimphy.h" +c + integer kmax2,imax2,kmin2,imin2 + real ftmax2,ftmin2 + integer kmax,imax,kmin,imin + real ftmax,ftmin +c + integer ncum + integer idcum(len) + integer len + integer nd + integer ndp1 + integer nl + integer minorig + integer nk1(len) + integer icb1(len) + real t1(len,nd) + real q1(len,nd) + real qs1(len,nd) + real u1(len,nd) + real v1(len,nd) + real gz1(len,nd) + real tv1(len,nd) + real tp1(len,nd) + real tvp1(len,nd) + real clw1(len,nd) + real h1(len,nd) + real lv1(len,nd) + real cpn1(len,nd) + real p1(len,nd) + real ph1(len,ndp1) + real ft1(len,nd) + real fq1(len,nd) + real fu1(len,nd) + real fv1(len,nd) + real tnk1(len) + real qnk1(len) + real gznk1(len) + real precip1(len) + real cbmf1(len) + real plcl1(len) + integer iflag1(len) + real delt + real cpd + real cpv + real cl + real rv + real rd + real lv0 + real g + real sigs ! SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE + real sigd ! SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT + real elcrit ! ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) + real tlcrit ! TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- +c CONVERSION THRESHOLD IS ASSUMED TO BE ZERO + real omtsnow ! OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW + real dtmax ! DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION +c A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC. + real damp + real alpha + real entp ! ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT FORMULATION + real coeffs ! COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION OF SNOW + real coeffr ! COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION OF RAIN + real omtrain ! OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN + real cu ! CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM TRANSPORT +c + real Ma(len,nd) +c +C +C *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) *** +C *** TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- *** +C *** CONVERSION THRESHOLD IS ASSUMED TO BE ZERO *** +C *** (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY *** +C *** BETWEEN 0 C AND TLCRIT) *** +C *** ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT *** +C *** FORMULATION *** +C *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** +C *** SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** +C *** OF CLOUD *** +C *** OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN *** +C *** OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW *** +C *** COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** +C *** OF RAIN *** +C *** COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** +C *** OF SNOW *** +C *** CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM *** +C *** TRANSPORT *** +C *** DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION *** +C *** A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC *** +C *** ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF *** +C *** APPROACH TO QUASI-EQUILIBRIUM *** +C *** (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) *** +C *** (DAMP MUST BE LESS THAN 1) *** +c +c Local arrays. +c + real work(ncum) + real t(ncum,klev) + real q(ncum,klev) + real qs(ncum,klev) + real u(ncum,klev) + real v(ncum,klev) + real gz(ncum,klev) + real h(ncum,klev) + real lv(ncum,klev) + real cpn(ncum,klev) + real p(ncum,klev) + real ph(ncum,klev) + real ft(ncum,klev) + real fq(ncum,klev) + real fu(ncum,klev) + real fv(ncum,klev) + real precip(ncum) + real cbmf(ncum) + real plcl(ncum) + real tnk(ncum) + real qnk(ncum) + real gznk(ncum) + real tv(ncum,klev) + real tp(ncum,klev) + real tvp(ncum,klev) + real clw(ncum,klev) +c real det(ncum,klev) + real dph(ncum,klev) +c real wd(ncum) +c real tprime(ncum) +c real qprime(ncum) + real ah0(ncum) + real ep(ncum,klev) + real sigp(ncum,klev) + integer nent(ncum,klev) + real water(ncum,klev) + real evap(ncum,klev) + real mp(ncum,klev) + real m(ncum,klev) + real qti + real wt(ncum,klev) + real hp(ncum,klev) + real lvcp(ncum,klev) + real elij(ncum,klev,klev) + real ment(ncum,klev,klev) + real sij(ncum,klev,klev) + real qent(ncum,klev,klev) + real uent(ncum,klev,klev) + real vent(ncum,klev,klev) + real qp(ncum,klev) + real up(ncum,klev) + real vp(ncum,klev) + real cape(ncum) + real capem(ncum) + real frac(ncum) + real dtpbl(ncum) + real tvpplcl(ncum) + real tvaplcl(ncum) + real dtmin(ncum) + real w3d(ncum,klev) + real am(ncum) + real ents(ncum) + real uav(ncum) + real vav(ncum) +c + integer iflag(ncum) + integer nk(ncum) + integer icb(ncum) + integer inb(ncum) + integer inb1(ncum) + integer jtt(ncum) +c + integer nn,i,k,n,icbmax,nlp,j + integer ij + integer nn2,nn3 + real clmcpv + real clmcpd + real cpdmcp + real cpvmcpd + real eps + real epsi + real epsim1 + real tg,qg,s,alv,tc,ahg,denom,es,rg,ginv,rowl + real delti + real tca,elacrit + real by,defrac +c real byp + real byp(ncum) + logical lcape(ncum) + real dbo + real bf2,anum,dei,altem,cwat,stemp + real alt,qp1,smid,sjmax,sjmin + real delp,delm + real awat,coeff,afac,revap,dhdp,fac,qstm,rat + real qsm,sigt,b6,c6 + real dpinv,cpinv + real fqold,ftold,fuold,fvold + real wdtrain(ncum),xxx + real bsum(ncum,klev) + real asij(ncum) + real smin(ncum) + real scrit(ncum) +c real amp1,ad + real amp1(ncum),ad(ncum) + logical lwork(ncum) + integer num1,num2 +c +c print*,'cpd en entree de convect2 ',cpd + nlp=nl+1 +c + rowl=1000.0 + ginv=1.0/g + delti=1.0/delt +c +c Define some thermodynamic variables. +c + clmcpv=cl-cpv + clmcpd=cl-cpd + cpdmcp=cpd-cpv + cpvmcpd=cpv-cpd + eps=rd/rv + epsi=1.0/eps + epsim1=epsi-1.0 +c +c Compress the fields. +c + do 110 k=1,nl+1 + nn=0 + do 100 i=1,len + if(iflag1(i).eq.0)then + nn=nn+1 + t(nn,k)=t1(i,k) + q(nn,k)=q1(i,k) + qs(nn,k)=qs1(i,k) + u(nn,k)=u1(i,k) + v(nn,k)=v1(i,k) + gz(nn,k)=gz1(i,k) + h(nn,k)=h1(i,k) + lv(nn,k)=lv1(i,k) + cpn(nn,k)=cpn1(i,k) + p(nn,k)=p1(i,k) + ph(nn,k)=ph1(i,k) + tv(nn,k)=tv1(i,k) + tp(nn,k)=tp1(i,k) + tvp(nn,k)=tvp1(i,k) + clw(nn,k)=clw1(i,k) + endif + 100 continue +c print*,'100 ncum,nn',ncum,nn + 110 continue + nn=0 + do 150 i=1,len + if(iflag1(i).eq.0)then + nn=nn+1 + cbmf(nn)=cbmf1(i) + plcl(nn)=plcl1(i) + tnk(nn)=tnk1(i) + qnk(nn)=qnk1(i) + gznk(nn)=gznk1(i) + nk(nn)=nk1(i) + icb(nn)=icb1(i) + iflag(nn)=iflag1(i) + endif + 150 continue +c print*,'150 ncum,nn',ncum,nn +c +c Initialize the tendencies, det, wd, tprime, qprime. +c + do 170 k=1,nl + do 160 i=1,ncum +c det(i,k)=0.0 + ft(i,k)=0.0 + fu(i,k)=0.0 + fv(i,k)=0.0 + fq(i,k)=0.0 + dph(i,k)=ph(i,k)-ph(i,k+1) + ep(i,k)=0.0 + sigp(i,k)=sigs + 160 continue + 170 continue + do 180 i=1,ncum +c wd(i)=0.0 +c tprime(i)=0.0 +c qprime(i)=0.0 + precip(i)=0.0 + ft(i,nl+1)=0.0 + fu(i,nl+1)=0.0 + fv(i,nl+1)=0.0 + fq(i,nl+1)=0.0 + 180 continue +c +c Compute icbmax. +c + icbmax=2 + do 230 i=1,ncum + icbmax=max(icbmax,icb(i)) + 230 continue +c +c +!===================================================================== +! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES +!===================================================================== +c +c --- The procedure is to solve the equation. +c cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. +c +c *** Calculate certain parcel quantities, including static energy *** +c +c + do 240 i=1,ncum + ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + & +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i) + 240 continue +c +c +c *** Find lifted parcel quantities above cloud base *** +c +c + do 300 k=minorig+1,nl + do 290 i=1,ncum + if(k.ge.(icb(i)+1))then + tg=t(i,k) + qg=qs(i,k) + alv=lv0-clmcpv*(t(i,k)-273.15) +c +c First iteration. +c + s=cpd+alv*alv*qg/(rv*t(i,k)*t(i,k)) + s=1./s + ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k) + tg=tg+s*(ah0(i)-ahg) + tg=max(tg,35.0) + tc=tg-273.15 + denom=243.5+tc + if(tc.ge.0.0)then + es=6.112*exp(17.67*tc/denom) + else + es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) + endif + qg=eps*es/(p(i,k)-es*(1.-eps)) +c +c Second iteration. +c + s=cpd+alv*alv*qg/(rv*t(i,k)*t(i,k)) + s=1./s + ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k) + tg=tg+s*(ah0(i)-ahg) + tg=max(tg,35.0) + tc=tg-273.15 + denom=243.5+tc + if(tc.ge.0.0)then + es=6.112*exp(17.67*tc/denom) + else + es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) + endif + qg=eps*es/(p(i,k)-es*(1.-eps)) +c + alv=lv0-clmcpv*(t(i,k)-273.15) +c print*,'cpd dans convect2 ',cpd +c print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd' +c print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd + tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg) + & /cpd +c if (.not.cpd.gt.1000.) then +c print*,'CPD=',cpd +c stop +c endif + clw(i,k)=qnk(i)-qg + clw(i,k)=max(0.0,clw(i,k)) + rg=qg/(1.-qnk(i)) + tvp(i,k)=tp(i,k)*(1.+rg*epsi) + endif + 290 continue + 300 continue +c +!===================================================================== +! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF +! --- PRECIPITATION FALLING OUTSIDE OF CLOUD +! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) +!===================================================================== +c + do 320 k=minorig+1,nl + do 310 i=1,ncum + if(k.ge.(nk(i)+1))then + tca=tp(i,k)-273.15 + if(tca.ge.0.0)then + elacrit=elcrit + else + elacrit=elcrit*(1.0-tca/tlcrit) + endif + elacrit=max(elacrit,0.0) + ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8) + ep(i,k)=max(ep(i,k),0.0 ) + ep(i,k)=min(ep(i,k),1.0 ) + sigp(i,k)=sigs + endif + 310 continue + 320 continue +c +!===================================================================== +! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL +! --- VIRTUAL TEMPERATURE +!===================================================================== +c + do 340 k=minorig+1,nl + do 330 i=1,ncum + if(k.ge.(icb(i)+1))then + tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k)) +c print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)' +c print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k) + endif + 330 continue + 340 continue + do 350 i=1,ncum + tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd + 350 continue +c +c +c===================================================================== +c --- NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS +c===================================================================== +c + do 360 i=1,ncum*nlp + nent(i,1)=0 + water(i,1)=0.0 + evap(i,1)=0.0 + mp(i,1)=0.0 + m(i,1)=0.0 + wt(i,1)=omtsnow + hp(i,1)=h(i,1) +c if(.not.cpn(i,1).gt.900.) then +c print*,'i,lv(i,1),cpn(i,1)' +c print*, i,lv(i,1),cpn(i,1) +c k=(i-1)/ncum+1 +c print*,'i,k',mod(i,ncum),k,' cpn',cpn(mod(i,ncum),k) +c stop +c endif + lvcp(i,1)=lv(i,1)/cpn(i,1) + 360 continue +c + do 380 i=1,ncum*nlp*nlp + elij(i,1,1)=0.0 + ment(i,1,1)=0.0 + sij(i,1,1)=0.0 + 380 continue +c + do 400 k=1,nlp + do 390 j=1,nlp + do 385 i=1,ncum + qent(i,k,j)=q(i,j) + uent(i,k,j)=u(i,j) + vent(i,k,j)=v(i,j) + 385 continue + 390 continue + 400 continue +c + do 420 i=1,ncum + qp(i,1)=q(i,1) + up(i,1)=u(i,1) + vp(i,1)=v(i,1) + 420 continue + do 440 k=2,nlp + do 430 i=1,ncum + qp(i,k)=q(i,k-1) + up(i,k)=u(i,k-1) + vp(i,k)=v(i,k-1) + 430 continue + 440 continue +c +c===================================================================== +c --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S +c --- HIGHEST LEVEL OF NEUTRAL BUOYANCY +c --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB) +c===================================================================== +c + do 510 i=1,ncum + cape(i)=0.0 + capem(i)=0.0 + inb(i)=icb(i)+1 + inb1(i)=inb(i) + 510 continue +c +c Originial Code +c +c do 530 k=minorig+1,nl-1 +c do 520 i=1,ncum +c if(k.ge.(icb(i)+1))then +c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) +c byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) +c cape(i)=cape(i)+by +c if(by.ge.0.0)inb1(i)=k+1 +c if(cape(i).gt.0.0)then +c inb(i)=k+1 +c capem(i)=cape(i) +c endif +c endif +c520 continue +c530 continue +c do 540 i=1,ncum +c byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl) +c cape(i)=capem(i)+byp +c defrac=capem(i)-cape(i) +c defrac=max(defrac,0.001) +c frac(i)=-cape(i)/defrac +c frac(i)=min(frac(i),1.0) +c frac(i)=max(frac(i),0.0) +c540 continue +c +c K Emanuel fix +c +c call zilch(byp,ncum) +c do 530 k=minorig+1,nl-1 +c do 520 i=1,ncum +c if(k.ge.(icb(i)+1))then +c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) +c cape(i)=cape(i)+by +c if(by.ge.0.0)inb1(i)=k+1 +c if(cape(i).gt.0.0)then +c inb(i)=k+1 +c capem(i)=cape(i) +c byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) +c endif +c endif +c520 continue +c530 continue +c do 540 i=1,ncum +c inb(i)=max(inb(i),inb1(i)) +c cape(i)=capem(i)+byp(i) +c defrac=capem(i)-cape(i) +c defrac=max(defrac,0.001) +c frac(i)=-cape(i)/defrac +c frac(i)=min(frac(i),1.0) +c frac(i)=max(frac(i),0.0) +c540 continue +c +c J Teixeira fix +c + call zilch(byp,ncum) + do 515 i=1,ncum + lcape(i)=.true. + 515 continue + do 530 k=minorig+1,nl-1 + do 520 i=1,ncum + if(cape(i).lt.0.0)lcape(i)=.false. + if((k.ge.(icb(i)+1)).and.lcape(i))then + by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) + byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) + cape(i)=cape(i)+by + if(by.ge.0.0)inb1(i)=k+1 + if(cape(i).gt.0.0)then + inb(i)=k+1 + capem(i)=cape(i) + endif + endif + 520 continue + 530 continue + do 540 i=1,ncum + cape(i)=capem(i)+byp(i) + defrac=capem(i)-cape(i) + defrac=max(defrac,0.001) + frac(i)=-cape(i)/defrac + frac(i)=min(frac(i),1.0) + frac(i)=max(frac(i),0.0) + 540 continue +c +c===================================================================== +c --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL +c===================================================================== +c + do 600 k=minorig+1,nl + do 590 i=1,ncum + if((k.ge.icb(i)).and.(k.le.inb(i)))then + hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k) + endif + 590 continue + 600 continue +c +c===================================================================== +c --- CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I), +c --- AT EACH MODEL LEVEL +c===================================================================== +c +c tvpplcl = parcel temperature lifted adiabatically from level +c icb-1 to the LCL. +c tvaplcl = virtual temperature at the LCL. +c + do 610 i=1,ncum + dtpbl(i)=0.0 + tvpplcl(i)=tvp(i,icb(i)-1) + & -rd*tvp(i,icb(i)-1)*(p(i,icb(i)-1)-plcl(i)) + & /(cpn(i,icb(i)-1)*p(i,icb(i)-1)) + tvaplcl(i)=tv(i,icb(i)) + & +(tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i,icb(i))) + & /(p(i,icb(i))-p(i,icb(i)+1)) + 610 continue +c +c------------------------------------------------------------------- +c --- Interpolate difference between lifted parcel and +c --- environmental temperatures to lifted condensation level +c------------------------------------------------------------------- +c +c dtpbl = average of tvp-tv in the PBL (k=nk to icb-1). +c + do 630 k=minorig,icbmax + do 620 i=1,ncum + if((k.ge.nk(i)).and.(k.le.(icb(i)-1)))then + dtpbl(i)=dtpbl(i)+(tvp(i,k)-tv(i,k))*dph(i,k) + endif + 620 continue + 630 continue + do 640 i=1,ncum + dtpbl(i)=dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i))) + dtmin(i)=tvpplcl(i)-tvaplcl(i)+dtmax+dtpbl(i) + 640 continue +c +c------------------------------------------------------------------- +c --- Adjust cloud base mass flux +c------------------------------------------------------------------- +c + do 650 i=1,ncum + work(i)=cbmf(i) + cbmf(i)=max(0.0,(1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i)) + if((work(i).eq.0.0).and.(cbmf(i).eq.0.0))then + iflag(i)=3 + endif + 650 continue +c +c------------------------------------------------------------------- +c --- Calculate rates of mixing, m(i) +c------------------------------------------------------------------- +c + call zilch(work,ncum) +c + do 670 j=minorig+1,nl + do 660 i=1,ncum + if((j.ge.(icb(i)+1)).and.(j.le.inb(i)))then + k=min(j,inb1(i)) + dbo=abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) + & +entp*0.04*(ph(i,k)-ph(i,k+1)) + work(i)=work(i)+dbo + m(i,j)=cbmf(i)*dbo + endif + 660 continue + 670 continue + do 690 k=minorig+1,nl + do 680 i=1,ncum + if((k.ge.(icb(i)+1)).and.(k.le.inb(i)))then + m(i,k)=m(i,k)/work(i) + endif + 680 continue + 690 continue +c +c +c===================================================================== +c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING +c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING +c --- FRACTION (sij) +c===================================================================== +c +c + do 750 i=minorig+1,nl + do 710 j=minorig+1,nl + do 700 ij=1,ncum + if((i.ge.(icb(ij)+1)).and.(j.ge.icb(ij)) + & .and.(i.le.inb(ij)).and.(j.le.inb(ij)))then + qti=qnk(ij)-ep(ij,i)*clw(ij,i) + bf2=1.+lv(ij,j)*lv(ij,j)*qs(ij,j) + & /(rv*t(ij,j)*t(ij,j)*cpd) + anum=h(ij,j)-hp(ij,i)+(cpv-cpd)*t(ij,j)*(qti-q(ij,j)) + denom=h(ij,i)-hp(ij,i)+(cpd-cpv)*(q(ij,i)-qti)*t(ij,j) + dei=denom + if(abs(dei).lt.0.01)dei=0.01 + sij(ij,i,j)=anum/dei + sij(ij,i,i)=1.0 + altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j) + altem=altem/bf2 + cwat=clw(ij,j)*(1.-ep(ij,j)) + stemp=sij(ij,i,j) + if((stemp.lt.0.0.or.stemp.gt.1.0.or. + 1 altem.gt.cwat).and.j.gt.i)then + anum=anum-lv(ij,j)*(qti-qs(ij,j)-cwat*bf2) + denom=denom+lv(ij,j)*(q(ij,i)-qti) + if(abs(denom).lt.0.01)denom=0.01 + sij(ij,i,j)=anum/denom + altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j) + altem=altem-(bf2-1.)*cwat + endif + if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then + qent(ij,i,j)=sij(ij,i,j)*q(ij,i) + & +(1.-sij(ij,i,j))*qti + uent(ij,i,j)=sij(ij,i,j)*u(ij,i) + & +(1.-sij(ij,i,j))*u(ij,nk(ij)) + vent(ij,i,j)=sij(ij,i,j)*v(ij,i) + & +(1.-sij(ij,i,j))*v(ij,nk(ij)) + elij(ij,i,j)=altem + elij(ij,i,j)=max(0.0,elij(ij,i,j)) + ment(ij,i,j)=m(ij,i)/(1.-sij(ij,i,j)) + nent(ij,i)=nent(ij,i)+1 + endif + sij(ij,i,j)=max(0.0,sij(ij,i,j)) + sij(ij,i,j)=min(1.0,sij(ij,i,j)) + endif + 700 continue + 710 continue +c +c *** If no air can entrain at level i assume that updraft detrains *** +c *** at that level and calculate detrained air flux and properties *** +c + do 740 ij=1,ncum + if((i.ge.(icb(ij)+1)).and.(i.le.inb(ij)) + & .and.(nent(ij,i).eq.0))then + ment(ij,i,i)=m(ij,i) + qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i) + uent(ij,i,i)=u(ij,nk(ij)) + vent(ij,i,i)=v(ij,nk(ij)) + elij(ij,i,i)=clw(ij,i) + sij(ij,i,i)=1.0 + endif + 740 continue + 750 continue +c + do 770 i=1,ncum + sij(i,inb(i),inb(i))=1.0 + 770 continue +c +c===================================================================== +c --- NORMALIZE ENTRAINED AIR MASS FLUXES +c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING +c===================================================================== +c +c + call zilch(bsum,ncum*nlp) + do 780 ij=1,ncum + lwork(ij)=.false. + 780 continue + do 789 i=minorig+1,nl +c + num1=0 + do 779 ij=1,ncum + if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))num1=num1+1 + 779 continue + if(num1.le.0)go to 789 +c + do 781 ij=1,ncum + if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))then + lwork(ij)=(nent(ij,i).ne.0) + qp1=q(ij,nk(ij))-ep(ij,i)*clw(ij,i) + anum=h(ij,i)-hp(ij,i)-lv(ij,i)*(qp1-qs(ij,i)) + denom=h(ij,i)-hp(ij,i)+lv(ij,i)*(q(ij,i)-qp1) + if(abs(denom).lt.0.01)denom=0.01 + scrit(ij)=anum/denom + alt=qp1-qs(ij,i)+scrit(ij)*(q(ij,i)-qp1) + if(scrit(ij).lt.0.0.or.alt.lt.0.0)scrit(ij)=1.0 + asij(ij)=0.0 + smin(ij)=1.0 + endif + 781 continue + do 783 j=minorig,nl +c + num2=0 + do 778 ij=1,ncum + if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) + & .and.(j.ge.icb(ij)).and.(j.le.inb(ij)) + & .and.lwork(ij))num2=num2+1 + 778 continue + if(num2.le.0)go to 783 +c + do 782 ij=1,ncum + if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) + & .and.(j.ge.icb(ij)).and.(j.le.inb(ij)).and.lwork(ij))then + if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then + if(j.gt.i)then + smid=min(sij(ij,i,j),scrit(ij)) + sjmax=smid + sjmin=smid + if(smid.lt.smin(ij) + & .and.sij(ij,i,j+1).lt.smid)then + smin(ij)=smid + sjmax=min(sij(ij,i,j+1),sij(ij,i,j),scrit(ij)) + sjmin=max(sij(ij,i,j-1),sij(ij,i,j)) + sjmin=min(sjmin,scrit(ij)) + endif + else + sjmax=max(sij(ij,i,j+1),scrit(ij)) + smid=max(sij(ij,i,j),scrit(ij)) + sjmin=0.0 + if(j.gt.1)sjmin=sij(ij,i,j-1) + sjmin=max(sjmin,scrit(ij)) + endif + delp=abs(sjmax-smid) + delm=abs(sjmin-smid) + asij(ij)=asij(ij)+(delp+delm) + & *(ph(ij,j)-ph(ij,j+1)) + ment(ij,i,j)=ment(ij,i,j)*(delp+delm) + & *(ph(ij,j)-ph(ij,j+1)) + endif + endif + 782 continue + 783 continue + do 784 ij=1,ncum + if((i.ge.icb(ij)+1).and.(i.le.inb(ij)).and.lwork(ij))then + asij(ij)=max(1.0e-21,asij(ij)) + asij(ij)=1.0/asij(ij) + bsum(ij,i)=0.0 + endif + 784 continue + do 786 j=minorig,nl+1 + do 785 ij=1,ncum + if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) + & .and.(j.ge.icb(ij)).and.(j.le.inb(ij)) + & .and.lwork(ij))then + ment(ij,i,j)=ment(ij,i,j)*asij(ij) + bsum(ij,i)=bsum(ij,i)+ment(ij,i,j) + endif + 785 continue + 786 continue + do 787 ij=1,ncum + if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) + & .and.(bsum(ij,i).lt.1.0e-18).and.lwork(ij))then + nent(ij,i)=0 + ment(ij,i,i)=m(ij,i) + qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i) + uent(ij,i,i)=u(ij,nk(ij)) + vent(ij,i,i)=v(ij,nk(ij)) + elij(ij,i,i)=clw(ij,i) + sij(ij,i,i)=1.0 + endif + 787 continue + 789 continue +c +c===================================================================== +c --- PRECIPITATING DOWNDRAFT CALCULATION +c===================================================================== +c +c *** Check whether ep(inb)=0, if so, skip precipitating *** +c *** downdraft calculation *** +c +c +c *** Integrate liquid water equation to find condensed water *** +c *** and condensed water flux *** +c +c + do 890 i=1,ncum + jtt(i)=2 + if(ep(i,inb(i)).le.0.0001)iflag(i)=2 + if(iflag(i).eq.0)then + lwork(i)=.true. + else + lwork(i)=.false. + endif + 890 continue +c +c *** Begin downdraft loop *** +c +c + call zilch(wdtrain,ncum) + do 899 i=nl+1,1,-1 +c + num1=0 + do 879 ij=1,ncum + if((i.le.inb(ij)).and.lwork(ij))num1=num1+1 + 879 continue + if(num1.le.0)go to 899 +c +c +c *** Calculate detrained precipitation *** +c + do 891 ij=1,ncum + if((i.le.inb(ij)).and.(lwork(ij)))then + wdtrain(ij)=g*ep(ij,i)*m(ij,i)*clw(ij,i) + endif + 891 continue +c + if(i.gt.1)then + do 893 j=1,i-1 + do 892 ij=1,ncum + if((i.le.inb(ij)).and.(lwork(ij)))then + awat=elij(ij,j,i)-(1.-ep(ij,i))*clw(ij,i) + awat=max(0.0,awat) + wdtrain(ij)=wdtrain(ij)+g*awat*ment(ij,j,i) + endif + 892 continue + 893 continue + endif +c +c *** Find rain water and evaporation using provisional *** +c *** estimates of qp(i)and qp(i-1) *** +c +c +c *** Value of terminal velocity and coeffecient of evaporation for snow *** +c + do 894 ij=1,ncum + if((i.le.inb(ij)).and.(lwork(ij)))then + coeff=coeffs + wt(ij,i)=omtsnow +c +c *** Value of terminal velocity and coeffecient of evaporation for rain *** +c + if(t(ij,i).gt.273.0)then + coeff=coeffr + wt(ij,i)=omtrain + endif + qsm=0.5*(q(ij,i)+qp(ij,i+1)) + afac=coeff*ph(ij,i)*(qs(ij,i)-qsm) + & /(1.0e4+2.0e3*ph(ij,i)*qs(ij,i)) + afac=max(afac,0.0) + sigt=sigp(ij,i) + sigt=max(0.0,sigt) + sigt=min(1.0,sigt) + b6=100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij,i) + c6=(water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij,i) + revap=0.5*(-b6+sqrt(b6*b6+4.*c6)) + evap(ij,i)=sigt*afac*revap + water(ij,i)=revap*revap +c +c *** Calculate precipitating downdraft mass flux under *** +c *** hydrostatic approximation *** +c + if(i.gt.1)then + dhdp=(h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i)) + dhdp=max(dhdp,10.0) + mp(ij,i)=100.*ginv*lv(ij,i)*sigd*evap(ij,i)/dhdp + mp(ij,i)=max(mp(ij,i),0.0) +c +c *** Add small amount of inertia to downdraft *** +c + fac=20.0/(ph(ij,i-1)-ph(ij,i)) + mp(ij,i)=(fac*mp(ij,i+1)+mp(ij,i))/(1.+fac) +c +c *** Force mp to decrease linearly to zero *** +c *** between about 950 mb and the surface *** +c + if(p(ij,i).gt.(0.949*p(ij,1)))then + jtt(ij)=max(jtt(ij),i) + mp(ij,i)=mp(ij,jtt(ij))*(p(ij,1)-p(ij,i)) + & /(p(ij,1)-p(ij,jtt(ij))) + endif + endif +c +c *** Find mixing ratio of precipitating downdraft *** +c + if(i.ne.inb(ij))then + if(i.eq.1)then + qstm=qs(ij,1) + else + qstm=qs(ij,i-1) + endif + if(mp(ij,i).gt.mp(ij,i+1))then + rat=mp(ij,i+1)/mp(ij,i) + qp(ij,i)=qp(ij,i+1)*rat+q(ij,i)*(1.0-rat)+100.*ginv* + & sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i)) + up(ij,i)=up(ij,i+1)*rat+u(ij,i)*(1.-rat) + vp(ij,i)=vp(ij,i+1)*rat+v(ij,i)*(1.-rat) + else + if(mp(ij,i+1).gt.0.0)then + qp(ij,i)=(gz(ij,i+1)-gz(ij,i) + & +qp(ij,i+1)*(lv(ij,i+1)+t(ij,i+1) + & *(cl-cpd))+cpd*(t(ij,i+1)-t(ij,i))) + & /(lv(ij,i)+t(ij,i)*(cl-cpd)) + up(ij,i)=up(ij,i+1) + vp(ij,i)=vp(ij,i+1) + endif + endif + qp(ij,i)=min(qp(ij,i),qstm) + qp(ij,i)=max(qp(ij,i),0.0) + endif + endif + 894 continue + 899 continue +c +c *** Calculate surface precipitation in mm/day *** +c + do 1190 i=1,ncum + if(iflag(i).le.1)then +cc precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000. +cc & /(rowl*g) +cc precip(i)=precip(i)*delt/86400. + precip(i) = wt(i,1)*sigd*water(i,1)*86400/g + endif + 1190 continue +c +c +c *** Calculate downdraft velocity scale and surface temperature and *** +c *** water vapor fluctuations *** +c +c wd=beta*abs(mp(icb))*0.01*rd*t(icb)/(sigd*p(icb)) +c qprime=0.5*(qp(1)-q(1)) +c tprime=lv0*qprime/cpd +c +c *** Calculate tendencies of lowest level potential temperature *** +c *** and mixing ratio *** +c + do 1200 i=1,ncum + work(i)=0.01/(ph(i,1)-ph(i,2)) + am(i)=0.0 + 1200 continue + do 1220 k=2,nl + do 1210 i=1,ncum + if((nk(i).eq.1).and.(k.le.inb(i)).and.(nk(i).eq.1))then + am(i)=am(i)+m(i,k) + endif + 1210 continue + 1220 continue + do 1240 i=1,ncum + if((g*work(i)*am(i)).ge.delti)iflag(i)=1 + ft(i,1)=ft(i,1)+g*work(i)*am(i)*(t(i,2)-t(i,1) + & +(gz(i,2)-gz(i,1))/cpn(i,1)) + ft(i,1)=ft(i,1)-lvcp(i,1)*sigd*evap(i,1) + ft(i,1)=ft(i,1)+sigd*wt(i,2)*(cl-cpd)*water(i,2)*(t(i,2) + & -t(i,1))*work(i)/cpn(i,1) + fq(i,1)=fq(i,1)+g*mp(i,2)*(qp(i,2)-q(i,1))* + & work(i)+sigd*evap(i,1) + fq(i,1)=fq(i,1)+g*am(i)*(q(i,2)-q(i,1))*work(i) + fu(i,1)=fu(i,1)+g*work(i)*(mp(i,2)*(up(i,2)-u(i,1)) + & +am(i)*(u(i,2)-u(i,1))) + fv(i,1)=fv(i,1)+g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1)) + & +am(i)*(v(i,2)-v(i,1))) + 1240 continue + do 1260 j=2,nl + do 1250 i=1,ncum + if(j.le.inb(i))then + fq(i,1)=fq(i,1) + & +g*work(i)*ment(i,j,1)*(qent(i,j,1)-q(i,1)) + fu(i,1)=fu(i,1) + & +g*work(i)*ment(i,j,1)*(uent(i,j,1)-u(i,1)) + fv(i,1)=fv(i,1) + & +g*work(i)*ment(i,j,1)*(vent(i,j,1)-v(i,1)) + endif + 1250 continue + 1260 continue +c +c *** Calculate tendencies of potential temperature and mixing ratio *** +c *** at levels above the lowest level *** +c +c *** First find the net saturated updraft and downdraft mass fluxes *** +c *** through each level *** +c + do 1500 i=2,nl+1 +c + num1=0 + do 1265 ij=1,ncum + if(i.le.inb(ij))num1=num1+1 + 1265 continue + if(num1.le.0)go to 1500 +c + call zilch(amp1,ncum) + call zilch(ad,ncum) +c + do 1280 k=i+1,nl+1 + do 1270 ij=1,ncum + if((i.ge.nk(ij)).and.(i.le.inb(ij)) + & .and.(k.le.(inb(ij)+1)))then + amp1(ij)=amp1(ij)+m(ij,k) + endif + 1270 continue + 1280 continue +c + do 1310 k=1,i + do 1300 j=i+1,nl+1 + do 1290 ij=1,ncum + if((j.le.(inb(ij)+1)).and.(i.le.inb(ij)))then + amp1(ij)=amp1(ij)+ment(ij,k,j) + endif + 1290 continue + 1300 continue + 1310 continue + do 1340 k=1,i-1 + do 1330 j=i,nl+1 + do 1320 ij=1,ncum + if((i.le.inb(ij)).and.(j.le.inb(ij)))then + ad(ij)=ad(ij)+ment(ij,j,k) + endif + 1320 continue + 1330 continue + 1340 continue +c + do 1350 ij=1,ncum + if(i.le.inb(ij))then + dpinv=0.01/(ph(ij,i)-ph(ij,i+1)) + cpinv=1.0/cpn(ij,i) +c + ft(ij,i)=ft(ij,i) + & +g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij,i) + & +(gz(ij,i+1)-gz(ij,i))*cpinv) + & -ad(ij)*(t(ij,i)-t(ij,i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) + & -sigd*lvcp(ij,i)*evap(ij,i) + ft(ij,i)=ft(ij,i)+g*dpinv*ment(ij,i,i)*(hp(ij,i)-h(ij,i)+ + & t(ij,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv + ft(ij,i)=ft(ij,i)+sigd*wt(ij,i+1)*(cl-cpd)*water(ij,i+1)* + & (t(ij,i+1)-t(ij,i))*dpinv*cpinv + fq(ij,i)=fq(ij,i)+g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij,i))- + & ad(ij)*(q(ij,i)-q(ij,i-1))) + fu(ij,i)=fu(ij,i)+g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij,i))- + & ad(ij)*(u(ij,i)-u(ij,i-1))) + fv(ij,i)=fv(ij,i)+g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij,i))- + & ad(ij)*(v(ij,i)-v(ij,i-1))) + endif + 1350 continue + do 1370 k=1,i-1 + do 1360 ij=1,ncum + if(i.le.inb(ij))then + awat=elij(ij,k,i)-(1.-ep(ij,i))*clw(ij,i) + awat=max(awat,0.0) + fq(ij,i)=fq(ij,i) + & +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-awat-q(ij,i)) + fu(ij,i)=fu(ij,i) + & +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i)) + fv(ij,i)=fv(ij,i) + & +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i)) + endif + 1360 continue + 1370 continue + do 1390 k=i,nl+1 + do 1380 ij=1,ncum + if((i.le.inb(ij)).and.(k.le.inb(ij)))then + fq(ij,i)=fq(ij,i) + & +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-q(ij,i)) + fu(ij,i)=fu(ij,i) + & +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i)) + fv(ij,i)=fv(ij,i) + & +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i)) + endif + 1380 continue + 1390 continue + do 1400 ij=1,ncum + if(i.le.inb(ij))then + fq(ij,i)=fq(ij,i) + & +sigd*evap(ij,i)+g*(mp(ij,i+1)* + & (qp(ij,i+1)-q(ij,i)) + & -mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv + fu(ij,i)=fu(ij,i) + & +g*(mp(ij,i+1)*(up(ij,i+1)-u(ij,i))-mp(ij,i)* + & (up(ij,i)-u(ij,i-1)))*dpinv + fv(ij,i)=fv(ij,i) + & +g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij,i))-mp(ij,i)* + & (vp(ij,i)-v(ij,i-1)))*dpinv + endif + 1400 continue + 1500 continue +c +c *** Adjust tendencies at top of convection layer to reflect *** +c *** actual position of the level zero cape *** +c + do 503 ij=1,ncum + fqold=fq(ij,inb(ij)) + fq(ij,inb(ij))=fq(ij,inb(ij))*(1.-frac(ij)) + fq(ij,inb(ij)-1)=fq(ij,inb(ij)-1) + & +frac(ij)*fqold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ + 1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*lv(ij,inb(ij)) + & /lv(ij,inb(ij)-1) + ftold=ft(ij,inb(ij)) + ft(ij,inb(ij))=ft(ij,inb(ij))*(1.-frac(ij)) + ft(ij,inb(ij)-1)=ft(ij,inb(ij)-1) + & +frac(ij)*ftold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ + 1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*cpn(ij,inb(ij)) + & /cpn(ij,inb(ij)-1) + fuold=fu(ij,inb(ij)) + fu(ij,inb(ij))=fu(ij,inb(ij))*(1.-frac(ij)) + fu(ij,inb(ij)-1)=fu(ij,inb(ij)-1) + & +frac(ij)*fuold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ + 1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij)))) + fvold=fv(ij,inb(ij)) + fv(ij,inb(ij))=fv(ij,inb(ij))*(1.-frac(ij)) + fv(ij,inb(ij)-1)=fv(ij,inb(ij)-1) + & +frac(ij)*fvold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ + 1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij)))) + 503 continue +c +c *** Very slightly adjust tendencies to force exact *** +c *** enthalpy, momentum and tracer conservation *** +c + do 682 ij=1,ncum + ents(ij)=0.0 + uav(ij)=0.0 + vav(ij)=0.0 + do 681 i=1,inb(ij) + ents(ij)=ents(ij) + & +(cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)-ph(ij,i+1)) + uav(ij)=uav(ij)+fu(ij,i)*(ph(ij,i)-ph(ij,i+1)) + vav(ij)=vav(ij)+fv(ij,i)*(ph(ij,i)-ph(ij,i+1)) + 681 continue + 682 continue + do 683 ij=1,ncum + ents(ij)=ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) + uav(ij)=uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) + vav(ij)=vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) + 683 continue + do 642 ij=1,ncum + do 641 i=1,inb(ij) + ft(ij,i)=ft(ij,i)-ents(ij)/cpn(ij,i) + fu(ij,i)=(1.-cu)*(fu(ij,i)-uav(ij)) + fv(ij,i)=(1.-cu)*(fv(ij,i)-vav(ij)) + 641 continue + 642 continue +c + do 1810 k=1,nl+1 + do 1800 i=1,ncum + if((q(i,k)+delt*fq(i,k)).lt.0.0)iflag(i)=10 + 1800 continue + 1810 continue +c +c + do 1900 i=1,ncum + if(iflag(i).gt.2)then + precip(i)=0.0 + cbmf(i)=0.0 + endif + 1900 continue + do 1920 k=1,nl + do 1910 i=1,ncum + if(iflag(i).gt.2)then + ft(i,k)=0.0 + fq(i,k)=0.0 + fu(i,k)=0.0 + fv(i,k)=0.0 + endif + 1910 continue + 1920 continue + do 2000 i=1,ncum + precip1(idcum(i))=precip(i) + cbmf1(idcum(i))=cbmf(i) + iflag1(idcum(i))=iflag(i) + 2000 continue + do 2020 k=1,nl + do 2010 i=1,ncum + ft1(idcum(i),k)=ft(i,k) + fq1(idcum(i),k)=fq(i,k) + fu1(idcum(i),k)=fu(i,k) + fv1(idcum(i),k)=fv(i,k) + 2010 continue + 2020 continue +c + DO k=1,nd + DO i=1,len + Ma(i,k) = 0. + ENDDO + ENDDO + DO k=nl,1,-1 + DO i=1,ncum + Ma(i,k) = Ma(i,k+1)+m(i,k) + ENDDO + ENDDO +c + return + end + Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/zilch.F =================================================================== --- LMDZ.3.3/branches/rel-LF/libf/phylmd/zilch.F (revision 265) +++ LMDZ.3.3/branches/rel-LF/libf/phylmd/zilch.F (revision 265) @@ -0,0 +1,13 @@ + subroutine zilch(x,m) +c +c Zero the real array x dimensioned m. +c + implicit none +c + integer m,i + real x(m) + do 1 i=1,m + x(i)= 0.0 + 1 continue + return + end Index: LMDZ.3.3/branches/rel-LF/physiq.def =================================================================== --- LMDZ.3.3/branches/rel-LF/physiq.def (revision 265) +++ LMDZ.3.3/branches/rel-LF/physiq.def (revision 265) @@ -0,0 +1,13 @@ +# +## $Header$ +# +# +# Automatically generated make config: don't edit +# + +OCEAN=force +VEGET=y +OK_journe=y +OK_mensuel=y +OK_instan=y +SCATTER=n